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1. INTRODUCTION 
A MOTIVATION FOR STUDY 


In recent years a novel approach to the field of digital signal processing has received 
significant attention. This new approach, the wavelet transform, has displayed the potential for 
applications in the fields of speech and image processing, in particular. Techniques based on the 


wavelet transform have been praised for, among other things, their efficiency of computation. 


Unfortunately for the electrical engineer specializing tn the field of digital signal 
processing, the concepts on which the wavelet transform have been based have been somewhat 
elusive. Pioneers of the field have been specialized in such diverse fields as geophysics, 
mathematical physics, and a branch of applied mathematics known as functional analysis. Not 
completely unrelated to the familiar Fourier transform methods, wavelet transforms extend the 
concept of signal decomposition through basis function expansion to a more general and abstract 
realm. Such abstract mathematical disciplines as set theory are invoked within the context of 


embedded vector spaces and their relationships to each other. 


Included among the purposes of this present study 1s an attempt to bridge the gap 
between the mathematician and engineer. It 1s intended to provide a rudimentary enough 
understanding of some of the concepts of functional analysis to facilitate a more in depth 
understanding. Additionally, a strony tie is demonstrated between wavelet-based methods of 


decomposition and areas within the field of electrical engineering which may be more familiar. 


Finally, an additional purpose of the present study 1s to evaluate wavelet transform 
methods for detection and identification applications. Wavelet processing methods will be 
applied to familiar signal processing problems such as the detection of signals in noise and the 
resolving of disparate frequency components of which a signal 1s composed. 

B. OUTLINE OF STUDY 

Chapter IT will introduce concepts from functional analysis which commonly appear in 
the literature on wavelet transforms. Not intended to provide an in depth introduction into the 
field of functional analysis, the chapter provides a few basic tools to support further study. In 
Chapter III, the basics of multirate system theory will be explored. In particular, the Quadrature 
Mirror Filter (QMF) bank for two and more channels will be demonstrated. Chapter IV will 
introduce the theory of multiresolution analysis. From the definitions introduced, Mallat's 
algorithm for the discrete wavelet transform will be introduced and shown to be equivalent to the 
QMEF bank. Two earlier multiresolution decompositions--the Laplacian pyramid and the a trous 
algorithm for the wavelet transform--will then be reviewed prior to extending the concepts of 
multiresolution analysis to general structures constructed from cascaded QMF banks. In the 
Chapter V, the basis functions of the wavelet transform and filters for QMF banks will be 
considered in greater depth. The chapter will include a discussion of two-scale difference 
equations and some basic filter design methods for QMF banks. Chapter VI will demonstrate for 
detection and identification applications the application of some multiresolution structures 


introduced in Chapter VI. Multiresolution structures will be compared with the periodogram 


decomposition from the perspective of computational efficiency, robustness with respect to 


noise, and the ability to resolve proximate spectral components. 


Il. RESULTS FROM FUNCTIONAL ANALYSIS 
A. INTRODUCTION 


In essence, wavelet transformation represents a recent development in the field of 
functional analysis. Consequently, a brief study of this branch of mathematics provides an 


appropriate starting point for a study of Wavelet Transform (WT) analysis. 


In mathematical terminology, a finctional is defined as [1] "a function that has a domain 
whose elements are functions, sets, or the like and that assumes [scalar] numerical values." 
Fourier transformation--an operation familiar to electrical engineers working in the field of 
signal processing--is an example of a functional analysis technique. The Fourter transform 
integral 

£(@)= <F { f(x) }(@) = f7 f(x) eF* dx (2.1) 
. a 


1s a functional which, through projection on a set of basis functions "e"’," maps a continuous 
function to a specific value for a given value of "w." Similarly, the Discrete Fourier Transform 
(DFT) is an example of a discrete functional. Fora discrete sequence x'=[x[0] x[1] ... x[M-1]], 


the DFT 


A=W OS ( 


tJ 
tJ 


r - -j2rkin 
where [W ].,, = Wy"" =e” mine 


maps the vector x to a discrete set of points X[k]. 


kan 
B. NORMS AND NORMED SPACES [2] 
Linear transformations generally involve the mapping of a function (or vector) from one 


vector space to another. The study of these transformations requires the capacity to quantify 


distances within these spaces. Measures for functional spaces possess analogies in 
geometry--such as length of a vector, distance between the points defined by two vectors, and 
the scalar product of two vectors. These measures, or norms, provide the ability to transform 
each element u in a vector space V into areal number. For a linear vector space V over the real 


number field R, the norm ||ul| for every element u € V satisfies the following conditions: 


|u|] = O, and ||ul| = 0 if and only if u = 0. (positivity) 
lox ul] = Jox|-|Jul] for ae R. (homogeneity) (3) 
Ju + vi] < |Jul] + || y|| (triangular inequality). 


Examples of norms which sometimes arise within the context of functional analysis 
follow. The most common norm is the L*-norm. Pythagorean theorem represents a special case 


of the L*-norm. For a function "x(t)" defined in the closed interval C[0, T], the 


L?-norm--denoted by ||x||,--is defined as 

IIxl], = i Ix(t)l? dt] (2.4) 

Another norm 1s the sup-norm or supremum norm. The supremum of a function or 

functional is defined as the least upper bound of a function. The concept of supremum is similar 
to the concept of maximum, however a function does not necessarily ever assume the value of 
the supremum. For example, consider a function f(x)=x?. If the domain of f(x) includes the 
closed interval [0, 1], then the extremum (and, incidentally, the supremum) of f(x) occurs for 
f(x)=1. If however, the domain of f(x) is restricted to the open interval (0, 1), then the function 


f(x) has no extremma. At its limit, f(x) approaches unity, however it can never assume that 


value. In this case, however, f(x) does possess a supremum, since f(x) is bounded above by 


unity. 


To applying the concept of supremum to define a sup-norm, consider the vector space 
V=C[0, t,], the set of all real-valued, once-differentiable, continuous functions of t in the closed 
interval [0, t,]. The sup-norm of x, ||x/|_ 1s defined as the supremum of the function "x(t)" on its 


domain, or 


xi, = sup{ |x(t)|: OSt Sty }. (2.5) 
The Lebesgue norm represents another norm which is related to the L’-norm. For a real 
number p € [0, 0) (1 R, the Lebesgue norm for the function u(t) defined on the interval [0, T] is 


defined as 


TT ‘ l/p. 
tall, = [fp huctyle dt 7 < 0, (2.6) 

If a norm can be defined in a given space, then that space can be characterized by that 

o An 
norm. for instance, if the Lebesgue norm is defined for a space of interest, than that space can 
be classified as a Lebesgue space. Furthermore, the space by which a norm 1s characterized 1s 
indicated in the subscript of the norm operator symbol. In some cases, such as for Lebesgue 
spaces, the subscript indicates the metric of the space. In others, such as sup-norm spaces, 
subscripts less indicative of the norm operation appear. In general, for some arbitrary, 
unspecified space U, the norm operator for that space is denoted by 

| + I. 
U 

This notation will be used throughout the remainder of this section to indicate general normed 


spaces. 


Inner product spaces are of significant interest to signal processing. An inner product 
space is defined to be a linear vector space on which an inner product can be defined. The 
concept of inner product is related to the definition of the norm of a space. For instance, the 
L’-inner product (u,v), is evaluated as 

(u, v)2 =|. u(x) > v"(x) dx (2m) 
for two vectors u, v € L7(Q). Unless otherwise indicated, throughout the remainder of this 
paper, inner products will be assumed to be L’ inner products and the corresponding 
distinguishing operator subscripts will be suppressed. Additionally, inner product operations 


possess the following properties: 


i (u, v) =(v, u) Symmetry 

pe (au, v)=Q(u, v) Homogeneity 

3 (u,tu,, v)=(u,, v) + (a, V) | Additivity (2.8) 
4, (u, u) > 0, and (u, u) = 0 if and only if u = 0. Positive Definite. 

5 I(u,v)| $ (uu) (y, v) = lull [lvl Cauchy-Schwartz inequality. 


Given an inner product space V, the concept of orthogonality 1s also important for 
characterizing the space. Two vectors u, v € V are said to be orthogonal if (u, v) = 0. 
Furthermore, it 1s possible to partition an inner product space into orthogonally complementary 
subspaces. For instance, consider the inner product space V,. If V,, 1s a subset of V, (denoted 
in mathematical symbology by V,, c V,), then the orthogonal complement W,., of V,_, is 


defined as 


W, = Ver = {ue V,: (uuv)=0 V ve Vf. (2.9) 
Additionally, the union of two orthgonally complementary subspaces V,, and W,_, to obtain a 


third is denoted V,, ® W,_, = V,. The concept of an orthogonally complementary subspaces 


proves critical when defining a multiresolution analysis. 


The Hilbert Space represents perhaps the most important inner product space for signal 


processing applications. An abstract Hilbert space is an inner product space which possesses the 
following characteristics [3]: 


l. Linearity--The operations of addition and of multiplication by real or complex numbers 


are defined for its elements: 


2 The metric of an Hilbert space is derived from its inner product. Consequently, for any 


two elements u and v, there 1s an associated real or complex number. 


Completenes2-If a sequence of elements {u,} satisfies the condition ||u, - u,,||/ 30 V 
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In 


m, n — 9, then there exists an element, u, such that |lu, - ul| ~ 0 Vn 7 ~. 


G LINEAR OPERATORS AND TRANSFORMATIONS [2] 


Many signal processing applications of interest involve linear transformations from one 
linear vector space U to another linear vector space V. If for some vector u in inner product 
space U, there 1s a corresponding vector v in image space V, then the operation T, by which u 
and v are related, constitutes the corresponding transformation. Mathematically, a mapping by T 


from space U to space V is denoted by T: U > V. 





With respect to some arbitrary linear transformation T: U — V, the following properties 


apply: 

le T{ou,t+Bu,} = aT{u,} + BT{u,} Linearity 

2a \|T{u,} - T{u,}|], <M lu, - u,||,, for some arbitrary M>0 Continuity 

3 \|T{u}||, <M |lull,, for some M > 0 Boundedness from Above (2.10) 
4. |T{u}|], 2 C |lul|,, for some C > 0 Boundedness from Below 


Two important consequences for linear operators result from the properties of (2.10). The first 
states that for some invertible, linear operator T: U — V which is bounded from below, 
IT" {v}I <[lvily/C. (2.11) 


Equation (2.11), proven in [2], expresses the Bounded Inverse Theorem. In other words, the 


inverse of a linear operation bounded from below is bounded from above. 


The second significant result stems from the linear property. Given T: U — V, for 


finite-dimensional vector spaces U and V, constituent vectors u and v, respectively, and bases 


{,, 9, ..-. O,} for U and {y,, W,, .... Ww} for V, then u and v can be represented as 
u= Si Ak DK 
i | gale) 
ae Bi Wk 


Now, writing down the form of the transformation, T: U - V, the terms in (2.12) are related by 


TW ™M3 


Biwis % on Tlox} = Thu}. (2.13) 


Since 9, € U and y, € V, the mapping of {,} — {w,} can be expressed as 


T {Ox} = th Wj. (2.14) 


In vector-matrix notation, {a,}, {B;} and {t,,} are related by 


B, {11 tio. ae QO, 

1 fo] 039 oes a 
Pi |e] ta 7 ak (2.15) 
hey tml tm2 Sue tmin Qn 


Equation (2.15) represents a linear transformation of the projections onto the bases of U and V. 


In signal processing, and other sciences involving representation in terms of linear 
transformations, mappings which are one-to-one are of significant interest. Given T: U > V, if 
for each ue U there exists a unique v € V, then the transformation T is defined to be 
isomorphic. Furthermore, for T{u} =v, given an isomorphic operator T, there exists a unique ~ 
inverse operator T”' such that u=T"'{v} [4]. Isormorphic operators are particularly useful 


because they can be inverted. With respect to signal processing, if a signal is decomposed by an 


invertible operator, then it can be reconstructed to recover the original signal or process. 


Before the final concept can be defined, elaboration on the definition of a functional is 
necessary. Linear functionals consist of a subset of linear transformations T: U - V in which 
the transform space V represents a scalar field [5]. In the case of functionals, the transform 
space V represents the dual or conjugate space of U. The dual space of U is assigned a special 
notation U. The notion of a dual space differs from that of a biorthogonal, or Riesz space 


which shall be addressed in a later section. 


Finally, the concept of adjoint operators occasionally plays a role in functional 


representation of sequences. Given a bounded, linear operator T: U — V for normed, linear 
spaces U and V, a linear functional } can be defined such that o(v)=o0(T{u})={(u). The resultant 
functional {--for u--is linear since ® and T are linear. It 1s possible, therefore, to define an 
adjoint operator T , such that &((u)={(T {v})=0(v). The functional f, by definition, lies in U’. 
Additionally, omitting the argument vectors, the transform T can be expressed as [=T’o. In the 
case of Hilbert spaces, an operator and its adjoint are related by inner products. Specifically, for 
u, vé i, a Hilbert space, 


(T{u}, v) =(v, T {v}). (2.16) 
D. REPRESENTATION IN INNER PRODUCT SPACES 


Representation theory consists of the theory of representing sequences or sets in terms of 
projections upon sets of vectors. The Riesz Representation Theorem provides the foundation for 
representation theory for functionals. If: represents a bounded linear operator in an Hilbert 
space if, the Riesz representation theorem [2] states that there exists a unique vector v, € #% such 
that 

(w)=(v, Ww) V ies HY 


This vector v, 1s called the representation for operator 0. 


The Fourier Series Theorem [2] provides concepts fundamental to representation of a function in 


a Hilbert space. For a countably infinite, orthonormal vector set {un}“, € #, then a series of the 


form £ &,u, converges if and only if & |E,|2 < ©. In this case, the series £ &, uy converges 
n=O n=0 n=0 


to the same element x irregardless of the ordering of the terms. Furthermore, the element x, the 


orthogonal set {u,} and the weighting factors {€} are related by 


G, = (x, u,) (2.17a) 
and 
x= D Equa. (2.17b) 
n=0 


By linearity, Parsevals Equality follows directly from (2.17): 

(y. x |= 5 ie Lie (y, us J (2.18) 

n=O 
Finally, if x = y, another relationship 1s obvious: 
7 < 2 
| x lJ =, x) = Z ex, ua)|?. (2.19) 
n=0 
At this point, consideration of additional concepts is justified. The concept of closure 

arises in representation theory [2]. A subset S of a normed space is said to be closed if it 
contains all its limit points. If S is a closed set, and if a sequence converges, then the limit to 
which the sequences converges 1s contained in S. The concept of denseness is related to closure. 


If S is an open set (not closed), the set of all additional points necessary to obtain a closed set 


including S is denoted by S. Consequently, the union of the sets SUS results in a closed set. A 
space S is said to be dense if for any vector v in S, there is an element in S which is arbitrarily 


close to v. Practically speaking, within a dense space S, it is, therefore, possible to define a 


representation with arbitrary precision if enough basis vectors are employed. 


Finally, in representation theory, the concept of completeness of an orthonormal set 


represents an important concept. Orthonormality, in the usual sense, refers to a set of vectors 
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{un }“., such that (u,, u,,) = 5,,, where 0,,, represents Kronecker's delta function. If an 


orthonormal basis font, € # is complete, then the only vector v € # such that (v, u,) = 0 V 


integer n € [1, N] is a vector such that ||v|| = 0. 


In representation by summing projections on a set of orthonormal basis functions U, one 


of three cases can exist [6]: 


le The set of basis functions 1s incomplete. 


to 


The set of basis functions is complete. 
S. Some subset U, of the basis functions constitutes a complete set where the 
complementary subset U, # ©. 


In case 1, it is not possible to obtain a complete representation from U. In other words, 
N 


n=l > 


(2.17) does not hold. Since U is incomplete; for the set of basis functions {u, } there 1s 
some non-zero vector v such that (v, u,) = 0 for any integer non [1, NJ]. The set of Rademacher 


functions [7] represents one such example. 


In case 2, a complete representation is obtainable and (2.17) does hold. Classical Fourier 


expansion represents an example for case 2. 


The third case is less common but does arise occasionally. In the third case, some 
arbitrary basis vector u, may lie in the closed linear span of all others in the set: 
lea Oy U,. f2rcu) 
ken 
Such an occurrence is referred to as a frame. In the case of a frame, (2.17) does not hold. Due 
to (2.20), evaluation of (2.17b) would result in a representation of vector x containing redundant 


information. Representation by Fourier expansion with overlapping subdomain basis functions 
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exemplifies the use of a frame. Specifically, representation by expansion on integer translates of 


polynomial splines is employed for some signal processing applications [9]. 


The concept of a frame originates from a theorem outlined by Riesz and Sz.-Nagy 
addressing biorthogonal systems [3]. The present term, however, was introduced by Duffin and 
Schaeffer in work related to nonharmonic Fourier analysis [8]. In contrast to expansion by sets 
of complete orthonormal bases, expansions by frames do not converge to a specific vector. 


Instead, such expansions converge to a specified range: 


All? <2], ¥)|" < Bik? Vke Z (2.21) 
where Z represents the set of all integers. The constant factors A, B € R are called the frame 
bounds. In order the operation to be invertible (or even, possibly, nontrivial), the requirement 
exists that A > 0 and B < ce. Daubechies [6] presents detailed descriptions of the procedures for 
calculating frame bounds ESE various Situations. 

Employing concepts from linear operator theory, the projection operation, T, of some 
vector, x onto a set of basis functions fo} VY me Zis characterized as a mapping from a 
Hilbert space # to the sequence of all square summable sequences (’(Z): 

T: H—> F(Z). 
Specifically, the projection operation T, defined to be the frame operator, is defined as 
(T{x}), = (O, x). (2.224) 


The adjoint frame operator, T , results in the expansion on the basis set {9,} 


es = 20k Ox. (2.22b) 


Combining (2.22a) and (2.22b), (2.21) can be expressed in its most abstract form as 


AlsTT<Bl 
where | denotes the identity operator. 


A decomposition by frames can be reconstructed exactly 1f a biorthogonal system, also 
occasionally referred to as a Riesz basis [9], is constructed. The form of such a reconstruction 


appears as 
X= 2 (Ok, X) Ok =Z (Ok, X) Ok. (2.23) 
In (2.23), the vector set {,} is referred to the biorthogonal basis of {,}. The basis function for 
the biorthogonal set is related to the fundamental set by 
je 1) oy (2.24) 
Occasionally, in literature, the Riesz basis will be referred to as the dual frame. This term 


should not be confused with the concept of dual or conjugate space previously discussed. 


Daubechies [6] describes and proves three generalities regarding biorthogonal systems. These 
include: 


1. The family LO. } uv: with O, = (T° T)7' 0, constitutes a frame with bounds B" and A’. 


2. The frame operator T associated with 0, is given by T =T(T*T)~! and satisfies 
ae ee 
and fees) 
T*T =T* T=1, 
or, equivalently 
CS 
be And finally, 


mr. =TT*. (2.26) 


IS 


HI. ASURVEY OF MULTIRATE SYSTEM THEORY 
A. INTRODUCTION 


Multiresolution analysis implementations frequently reduce to structures composed of 
multirate system building blocks. Consequently, prior to addressing the construction of 
multiresolution structures, multirate system theory will be briefly reviewed. Section B will 
address basic multirate system operations and representations. The simplest multirate 
system--the two-channel Quadrature Mirror Filter (QMF) bank--will be discussed in Section C. 
In Section D, the results from the preceding sections will be extended to an arbitrary number of 
channels. 

B. BASIC MULTIRATE SYSTEM OPERATIONS 

Multirate systems are comprised of three fundamental operators [10]--decimators, 

expanders, and linear (usually Finite ee Response (FIR)), digital filters. In this section, 


time- and Fourler-domain consequences of decimation and expansion will be demonstrated. 


Decimation consists of subsampling a discrete sequence, retaining only samples at 
integer intervals. Figure 3.1 presents the block diagram symbol for a decimator. In 


mathematical 


1M 


Figure 3.1--Block diagram for M-fold decimator. 


notation, decimation of a discrete sequence by a factor of "M" 1s denoted by 


X al )=x( Mr). 
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To illustrate the consequences of decimation, consider the bandpass, transient sequence 
plotted below in Figure 3.2. Figure 3.2 displays a lineplot of a 128-point sequence constructed 


from the real part of an exponentially-modulated, Kaiser windowed, sinc (sine-over-argument) 





pulse: 
e ane ; 
sin(nn’8 | ot Le eal | 
Re{s(n)} = —==— ——————  Re{e!™™3} (3.2) 
lo{ 7? ) 
8.15 
Q.1 
® 
=, 
Te 
> 
wy 
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% 
Ww 
"+65 45 -25 § 15 35 55 
Sample sequence number, n 
Figure 3.2--Time plot of discrete sequence described by equation (3.2). 
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Figure 3.3--Plot of sequence obtained from decimation by a factor of two of sequence (3.2). 
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Decimation by a factor of two, s,,(n), of the sequence in Figure 3.2, produces sequence 
plotted in Figure 3.3. The length of sequence s,,(n) is one half the length of sequence s(n). 


Furthermore, the sequences are related by 


Sioa), (335) 


To consider the Fourier-domain consequences of the decimation operation, the Z-transform of 


the decimated sequence must first be evaluated. For the decimated sequence, 


s(2n)z™ 


S(mizee* 


S 12(Z) 


* 4 
5 (3. 4) 


Evaluation of the second summation in (3.4a) requires the definition of a sequence 1),(n) whose 
values are unity for even elements and zero for odd elements [11]: 

n,(n) = [1 + (-1)"]/2 (3.5) 
Inserting (3.5) into (3.4) produces 


Sio(z) =P siz 


+E s(n) (1 +(-1)" 20? | (3.6) 
=| S(z'?) a S(-z"/?) ] 


The effect of decimation on the spectral content of the sequence s(n) 1s obtained from evaluation 


of S,,(z) on the unit circle: 


Si2(ei®) = L[S(ei?)+S(-e) 9?) | 


, | : 3.7 
PS (e yen o(e! (2m) | ( ) 


Decimation, therefore, causes the frequency spectrum of a frequency to become spread become 


by a factor of two. Additionally, the location of a spectral peak or any other distinguishable 
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feature will be translated by a factor of two with respect to its original spectral location. The 


Fourier transform indicated by (3.7) is a 4-%-periodic member of L’([0, 4-7]). 


— Original bandpass sequence 
—ODecinated bandpass sequence 
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Figure 3.4--Plots of magnitudes of S(e®) and S$ ie ae 


Figure 3.4 compares the Fourier transforms of a sequence S(e!®) and its decimated 


version S,,(e’°). The plots in Figure 3.4 illustrate the test sequence power distribution relative 
to its sampling frequency. The changes to sequence power distribution which occur from 
decimation are a consequence of a change in the equivalent sampling frequency. In discarding 


each separate sample, as is done in the decimation operation, the resultant sequence is equivalent 


to sampling the original sequence at a sample frequency of f,,, = 0.5-f,. 


As predicted by (3.7), decimation causes a spreading of a bandpass about its center 
frequency. Furthermore, distance between the center frequency f.,, of the decimated sequence 
and w=0 has been doubled with respect to that of the original sequence. In the example 


illustrated, the original sequence s(n) was constructed to be an analytic, bandpass process whose 
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center frequency f. was related to its sampling frequency f, by f, = f£/6 and whose bandwidth 
B=0.125-f.. The decimation operation shifts the center frequency f.,, = f/3 and enlarges the 
bandwidth to B,, = .25 f, or by a factor of two. The spreading and shifting effects occur because 
of the subsampling nature of the decimation operation: Decimation results in a change in the 
sampling frequency. The sequence s,,(n) is equivalent to s(n) sampled at one half of its original 


sampling frequency. The peak magnitude of S,,(e’®) is one half the magnitude of S(e”) because 


of the factor of one half introduced by the sifting sequence (3.5). 


—— 


Figure 3.5--Block diagram of an expansion operator. 
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Figure 3.6--Discrete plot of expanded version of s,,(n). 


The operation of expansion (or upsampling) represents, in a sense, the inversion of the 
decimation operation. Figure 3.5 presents the signal processing block diagram symbol for an 


expander. The mathematical notation for a two-fold expansion operation 1s given by: 
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s(n/2) V even ’n’ 
St2(n) me) 


0 vy odd Tl (3.8) 
£(1+(-1)* ) stn) 


Consequently, expansion entails the insertion of a zero between each sample in a sequence. 
Figure 3.6 presents a discrete line plot of the expansion of the sequence s,,(n). Expansion 
restores the length of the decimated sequence to that of the original sequence. However, all 


odd-number elements of the sequence s,,(n) are zero. 


Considering the Z-transform and Fourier transform of the expanded sequence provides 
further insight into the relationship between the expansion and decimation operations. The 


Z-transform of s;,(n) is evaluated as: 


St2(z) = $E(1+(-0"] s(2)2 
| (1 +(-1)** sin) | 
S 


= = 
: (8:9) 
= x (n) 
= s(z 
Consequently, in the Z-domain, expansion of the decimated sequence s,,(n) produces 
[Si2]42(z) = =[S(z) + S(-z)]. (3.10a) 


In other words, expansion of a decimated sequence reproduces the original sequence plus an 
additional term. This term is referred to as the aliasing term [12]. Evaluating (3.10a) on the 


unit circle yields: 


Za 


[Sir]n(e* | =[ S(e!®) + S(-e!®) ] 


| | (3.10b) 
= Si (6!2) Sie”) | 


a.6 —Decimated sequence 
} -—Expanded, decimated sequence 
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Figure 3.7--Magnitude plots of S,,(e'®) and [S,,],(e’”). 


The aliasing term, therefore, consists of the original term shifted by 7 radians. The magnitude 


plots of S,,(e°) and [S,,],,(2°) presented in Figure 3.7 illustrate this occurrence. The left side 


of the plot illustrates the restoration of the decimated sequence to its original center frequency 
and bandwidth while the nght side of the plot demonstrates the generation of an aliasing term. 
The aliasing term illustrated in Figure 3.7 can be reduced through the use of a linear, lowpass 


filter. Such a filter is often referred to as an interpolation filter. 


Figure 3.8 illustrates the time-domain results of applying an interpolation filter to the 
sequence plotted in Figure 3.6. The sequence in Figure 3.8 differs from the original sequence by 
a normalized mean-squared error of -11.52 dB where the normalized mean-square error 1s 


defined as 
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Net 42 
SD (sin)-<n) | 


is | s(n)—S(n) |” n=0 


2 N-] 
| s(n) | 2 s?(n) 


n=0 


(3.11) 





The interpolation filter whose response is shown in Figure 3.9 was constructed using a procedure 
outlined by Vaidyanathan [12]. A zero-phase, half-band, FIR filter was constructed using the 
McClellan-Parks algorithm and its characteristic polynomial factored. The minimum-phase 
zeros of the original polynomial were expanded to form an new characteristic polynomial for the 


filter whose frequency response appears in Figure 3.9. 
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Figure 3.8--Superimposed plots of original sequence Re{s(n)} and of real part of interpolated, expanded, 


decimated sequence, Re{h,*[s 1,0) } where h(n) represents the impulse response of the interpolation 
filter. 


The source of the error between the original sequence and its interpolated versions 1s 
readily apparent in Figure 3.9. The spectral content of the interpolated sequence almost exactly 
coincides with that of the original sequence in the region below the Nyquist frequency. 
However, since the interpolation filter consisted strictly of real coefficients, its frequency 


response included an image in the half of the spectrum above the Nyquist frequency. The 


location of the image coincided with the location of the aliasing term generated by expansion. 
As aresult of partial transmission of the aliasing component, the interpolated sequence whose 
spectrum is plotted in Figure 3.9 is no longer an analytical sequence. This accounts for the loss 


of symmetry evident in the plot of the interpolated sequence presented in Figure 3.8. 


St © 8 
ry” —Original sequence 
— Interpolated sequence 
+> Interpolatar filter response 


* Seay 
Va \ 


Normalized Magnitude (dB) 





Q a. @.2 8.3 Q.4 8.5 8.6 6.7 a.3 2.9 l 
Digital frequency, (multiple of f ) 
5 


Figure 3.9--Superimposed plots of spectral content of original sequence and interpolated sequence and of 
interpolator frequency response. 


C. TWO-CHANNEL QUADRATURE MIRROR FILTER BANKS 


The primary objective of the present section 1s to introduce the Quadrature Mirror Filter 
(QMF) bank. QMF bank theory 1s a well-developed area within the field of electrical 
engineering and can be understood with a grasp of basic linear systems theory. The purpose for 
the study of QMF banks in this chapter is to prepare for the construction of multiresolution 
structures in the next chapter. In Chapter IV, Mallat's algorithm for the discrete wavelet 
transform will be shown to be exactly equivalent to the QMF bank. Thereafter, multiresolution 


analysis structures will be constructed from QMF banks of various numbers of channels. 
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Figure 3.10--Block diagram of two-channel, Quadrature Murror Filter (QMF) bank. (After [12]) 


Having considered the basic building blocks of a linear, multireate system, attention in 
this section will be directed towards assembling those components into a basic system. 
Two-channel, quadrature mirror filter (QMF) banks represent the most basic structures for 
transmultiplexers, sub-band coders, and discrete wavelet transforms. Figure 3.10 presents a 
block diagram of a basic, two-channel QMF bank. The vertical, dotted line in Figure 3.10 
divides the system into analvsis and synthesis banks. Each of the filters in the structure is a 
half-band filter. The analysis section divides the signal, by frequency, into two channels. The 


synthesis bank then recombines the channels and generates an approximation x(n) of the original 


signal. In general, x(n) will differ from the original sequence because of three sources of error 
[13]: aliasing, amplitude distortion, and phase distortion. The Z-transform analysis presented in 


the previous section permits precise characterization of these errors. 


If a filter bank is designed such that each of these errors is exactly cancelled, then the 
structure 1s called a perfect reconstruction filter bank. More specifically, the perfect 
reconstruction criterion 1s expressed mathematically as 


x(n) =c- x(n—no) (3.12) 
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for some non-zero constant c and positive integer n,. In other words, for a perfect reconstruction 
filter bank, the reconstructed sequence will differ from the original sequence only by a constant 


factor and a delay. 


Beginning Z-transform analysis of the structure illustrated in Figure 3.11 [12], the output 
of the synthesis filter for channel k (k = 0, 1) 1s expressed as 
X,(Z)=H,(z) X(Z). (Sale 


By (3.6), the decimator output is described by 


Viz) = 5[Xx(2"?) + Xi(-z!?)] (3.14) 
= ALHy(2'?) X(2!) + Hy(-2") X(-27")] | 
Applying (3.9) and (3.10), 
Ux(z) = Vi(z?) (3.15) 


= +[Hy(z) X(z) +Hx(-z) X(-z)] 
describes the decimator output. Finally, the reconstruction component for channel k is given by 


X(z) = F(z) Ux(z) 


= +F,(z)[Hx(z) X(z)+Hi(-z) X(-z)] 


(3.16) 


Since the reconstructed signal is simply the sum of the synthesis filter bank channels, the overall 


transfer function for the filter bank, in matrix form, becomes 


; Ho(zie Elaitz) F(z) 
X(z)=+ X(z) X(-—z : : ae 
= aL ee@ ] BGea are 2) | F ,(Z) | — 


Now, (3.17) can be expanded and expressed as 


X(z) =+ F(z) Ho(z) + F(z) Hi(z) | X(z) + = | Fo(z) Ho(—z) + F(z) Hi(-z) | X(-z). 
2 2 
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When written in this manner, two categories of components are evident: the desired 
reconstructed components composed of X(z) and its factor, and an aliasing term consisting of 


X(-z) and its factor. Evaluated on the unit circle, the spectrum of the aliasing terms consist of a 
replica of the spectrum of the original sequence only shifted by z. 
In order to separately represent the components in the QMF bank output in terms of the 


reconstructed signal component and an undesired aliasing component, it is possible to represent 


(3.17) as a double-input, single output transfer function: 


T(z) |_a} Holz) Hi(z) Fo(Z) |_a}  Fo(z) Ho(z) + Fi(z) Hi(z) (3.18) 
A(z) | °| Ho(-z) Hi(-z) }| Fi(z) | *| Fo(z)Ho(-z)+Fi(z)Hi(-z) | 
To satisfy the perfect reconstruction criteria, it 1s necessary that 
A(e“)=0 Vi @ (3.19a) 
and 
Tei?) = eFom| Tigiey | =c-erion (3.19b) 





where c constitutes a positive constant and n, Is a positive integer. 


N 
Vaidyanathan [12] points out that, given a filter Ho(z) = 2 ho(n)z™, the remaining 


filters required to construct a two-channel, perfect reconstruction QMF bank are obtained from 


H(z) = -zHj(-1/z*) 
= -Nu*(|/>" 
Roz) = z  Hi(liz) (3.20a) 
and 
Eee = oz Hi (i/z") 
coupled with the constraint that 
: 2 | 2 
| Ho(ei®) | +| Hole) |” =c (3.20b) 
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for some positive constant c (not-necessarily the same constant as in (3.19b)) satisfy the perfect 
reconstruction criterion. A set of filters characterized by (3.20b) is said to be power 
complementary. Given (3.20), the problem of designing a two-channel, perfect reconstruction 
QMF filter bank reduces to the problem of selecting a valid analysis filter h(n). Some of the 


details of selecting good filters for these applications will be addressed in Chapter V. 


At this juncture, it is appropriate to introduce some additional concepts. The matrix 





H(z) = Ho(z) H(z) (3.21) 
Ho(-z) Hi(-z) 
constitutes the alias compensation matrix. If 
H#(e)®) H(e!®) =d I, (3.22) 


where the superscript ® denotes the transpose-conjugate of a matrix, I represents an identity 
matrix, and d is a positive constant, then the matrix H(e”) is said to be unitary. By definition, 
unitary operators are operators which, when applied to their inverses, produce an identity 
operation. Within the context of linear algebra, a matrix 1s unitary if each of its columns 1s 
linearly independent from all of the other columns in the matrix. Consequently, for a QMF bank 
with an alias compensation matrix which is unitary, the aliasing component of the output 1s 
linearly independent of the reconstructed component and the two components are, therefore, 
separable. 

Using the ti/de notation introduced by Vaidyanathan [12] where H(z)=H4( 1/z"), and 
where each element of H(z) is stable, if 


H(z) H(z) = d] (3235 
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then H(z) is said to satisfy the paraunitary property and the system characterized by H(z) can be 
described as a paraunitary system.. A paraunitary Z-transform matrix will be unitary when 
evaluated on the unit circle. Furthermore, if the coefficients of H(z) are all real, then H(z) is 
lossless bounded real. Additionally, Vaidyanathan lists three important properties for 


paraunitary systems: 


i The determinant of a square, FIR, paraunitary system produces an allpass polynomial. 


That 1s, 
det{H(z)} =az*,K 20,a¥#0.. (3.24) 
2. Paraunitary systems are power complementary. If h(?®)=[H,(e°) H,(e°)]? 


IH,(e)? + |H,(e*)P = h(el®) h(el®) =c Vw (3.25) 





>. The submatrices of H(z) are paraunitary. 


Another concept which becomes useful in multirate system theory is the polyphase 
representation. This representation will be introduced for a two-channel system here and 
extended in the next section to a structure entailing an arbitrary number of channels. It is 


possible to express the Z-transform of a sequence h,(n) as 


Hy(z) = Xhy(n)z™ 
‘ (3.24) 
= Zhy(2n)z~** + z7! Zhy(2n+ 1)z728 
If subsets of the sequence h,(n) are defined as 
exo(n) = hy(2n) | (3.25) 
exi(n) = hy(2n+1) 
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then it is possible to represent the two-channel, multirate system described in this channel in an 


additional manner. A z-dependent column vector h(z) can defined such that 


_| Ho(z) |_| Eoo(z*) Eos(z?) ee ee | 
n= | H(z) |-| E ico Zia ee zene) I i |-80 | Zia | — 


where 


bah Eoo(Z*) Eo3(z?) | 


Ei ote nee) 
Furthermore, upon comparing (3.26a) and (3.21), the following relationship arises: 


H7(z)=| h(z) h(-z) |. (3.26b) 


Therefore, in terms of the polyphase matrix E(z’) of (3.26a), (3.26b) becomes 


l l Lao fa) 
; h(z) h(—z) I-| a Jeed-| — I va fee. 


Consequently, a simple relationship between H(z) and E(z) exists: 


We =Ee)| - | 7 ‘ | Gm 


Obviously, the right-most matrix on the right-hand side of (3.27) is simply the transformation 
W, matnx for a two-point DFT. The middle matrix is simply a diagonal matmx of delays and 


can be expressed as 


] aad 
D3@)| 0 ool } 


Finally, it is not difficult to demonstrate that if the alias compensation matrix 1s 


paraunitary, then so is the polyphase matrix: 
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or equivalently, 


H(z) H(z) = E((z~)*) D2(1/z) W2 W# D2(z) E7(z?) 
= Pee) EB (z-) ' (3.28a) 
= d I 
E(z2) E(z2) =d I. (3.28b) 


Equation (3.28a) follows since W, Ww," = Land D,(z) D,(1/z) = 1. Equation (3.28b) follows 


from (3.28a) by making the substitution € =z ' followed by evaluating the complex conjugate of 


each side of the equation. 


Vaidyanathan delineates four important properties of paraunitary filter banks satisfying 


ieez0a) [12]: 


tI 


Lad 


Where N represents the order of the filters from which the structure 1s constructed, filter 


banks satisfying (3.20a) result in perfect reconstruction with 
| x(n) =0.5x(n—N); (3.29a) 
The analysis filters are power complementary and satisfy the condition 
| Hi(ei®) [=| Hotei) |; (3.29b) 
The synthesis filters are also complementary and, furthermore, they also satisfy 
| Fe(ei®) J =| Axle?) |; (3.29¢) 


and, all filters have order N=2J+1 for integer J. This condition ensures even filter 


lengths. 


To demonstrate the operation of a two-channel, QMF bank, a simple, 256-point sequence 


was constructed. The sequence was comprised of two sinusoids windowed by a complicated, 
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Figure 3.1 1--Time-domain plot of 256-point test sequence generated using (3.30). 
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Figure 3.]12--Plot of normalized power spectral density of 256-point test sequence generated using (3.30). 


exponential window: 
s(n) “(1 — e764 PPA gel \ Sim G n | += COS (2s n | (3.30) 


A time plot of the 256-point test sequence (3.30) is presented in Figure 3.11 and the power 


spectral density of the test sequence appears in figure 3.12. The test sequence contains harmonic 


a2 


components at digital frequencies of 7:7/20 (= 0.175-f,) and 17-2/20 (= 0.425-f.). The 


components were selected such that one appears in each half of the frequency spectum. 
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Figure 3.13--Impulse responses for 57°"-order FIR filters used in two-channel QMF analysis bank. 


Figures 3.13, 3.14 and 3.15 characterize the filter bank to which test sequence (3.30) was 
applied. The impulse response of the analysis filters for each channel are plotted in Figures 
3.13a and 3.13b. The low-frequency channel analysis filter, H,(z), a 57"-order FIR filter, is a 
minimum-phase spectral factor of a zero-phase, lowpass, half-band filter designed using the 
McClellan-Parks technique. From H,(z), the low-frequency channel synthesis filter F(z) and 


both high-frequency channel filters were designed using (3.20a). 


The equivalent transfer function, t(n), whose Z-transform was defined in (3.18), is 
plotted in Figure 3.14. It has been asserted in (3.12) and (3.19b) that, for perfect reconstruction, 
t(n) will be of the form 


t(n)=5(n-N), (3.31) 
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where N is the order of the filters from which the QMF bank was constructed. Superimposed 
with a dashed-line plot, is a plot of the sequence t(n)-8(n-N). The axis for the superimposed plot 
is graduated on the nght-hand side of the plot. From the supemmposed plot, it is evident that the 


sequence t(n) approaches the form of (3.31) with a maximum error of less than 2*10°. 
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Figure 3.14--Superimposed plots of equivalent impulse response for QMF bank, t(n) (defined as in (3.18)), and of 


t(n)-d(n-N,_.,). 
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Figure 3.15--Superimposed plots of frequency responses (in dB) of low- and high-frequency filters for QMF bank. 
Also superimposed is a plot of |H,(e”)|’+{H,(e”)/’. 


34 





= 
15 Magia 7 
e ; ? is 
? Pie ilig.s 
—_ A 4 ‘ tos ee 
Se 12 ree i= ei! ig 
— iter = 6 <i e ees’) 
un tae: *: ferme, 
: z ae | ees 
Ue ai i: Ac) sae ? 
q Staats ie eins 2S 
3 es eo ie tS 
Ty eran ee ee tet ee OF Fm» 
—a a es ey ee i eo Ss lees 10 @ ‘a! a aie ee me oe NL nm me 
ra 2 EN - CRE ects cet eet! st ue 8! eraser atm Na - - 1% eee: 
> * ee MG OC Oe ny Oy Ml ng 0 Te ahians oe" v4e 
ape a © Oc eee go SR Tay, es een ro ene} e oe 4 
‘ Lee ee Ye eae , ' hoe ' . . ‘ oe ty e 
qu 
— 
eS. 
= 
%~ 
Ww 





g 20 42 63 89 122 128 148 169 
Sample sequence number, n 


Figure 3.16--Time-domain plot of low-frequency channel decimator output for 256-point test sequence. 
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Figure 3.17--Time-domain plot of high-frequency channel decimator output for 256-point test sequence generated 
by (3.30). 


The frequency responses of the filters are plotted in Figure 3.15. Only two frequency 
responses are shown because of the condition described in (3.29c). Also superimposed on 
Figure 3.15 is a plot of |H,(e)|"+/H,(e)/?. This last plot demonstrates the degree to which the 


filter bank satisfies the power complementary property, (3.20b). 


BD 


Time plots of the low- and high-frequency channel decimator outputs are presented in 
Figures 3.16 and 3.17. As a consequence of the decimation-induced translation of the spectral 
peaks for each component, the plots appear to reflect harmonics whose frequencies are quite 
close together. Additionally, due to the shape of the filter, the flat segment in the trailing edge 
of the sequence plotted in Figure 3.16 converges sharply to zero. The flat segment in Figure 
3.17 is on the leading edge of the sequence. Since the synthesis filter for channel is the time 
reversal of the analysis filter, the relative locations of the flat regions will be reversed, resulting 
in a flat segment on either side. This represents the source of the delay. The difference occurs 
since, by (3.20a), the low-frequency and high-frequency channel analysis filters are modulated 
time-reversals of each other. Also worth ice the relative magnitudes of the low-frequency 
and high-frequency channel decimator outputs have maintained the proportions established in 
(3.30). As expected, the decimator output of the low-frequency in figure 3.16 channel maintains 
a peak amplitude greater, by a factor of four, than the peak if the envelope for the 


high-frequency channel in figure 3.17. 


Figure 3.18 presents superimposed, normalized power spectral density plots for the 


decimator outputs and the original sequence. As discussed previously, decimation changes the 
sampling frequency. The sampling frequency f,,, for a decimated sequence 1s related to the 


sampling frequency f, for the original frequency by 
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As expected from (3.7), the peak at 7/20 (= 0.175 f.) has been translated to 72/10 (= 0.35 Loe 
The peak which was originally at 17x/20 (= 0.425 f,) was shifted to 17m/10. Due to aliasing, the 


spectral peak located at 177/10 appears at 207/10 - 172/10 = 37/10 (= 0.15 f,),). Additionally, 


each of the peaks has been widened. 
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Figure 3.18--Supenmposed plots of normalized power spectral densities of the high- and low-frequency channel 
decimator outputs and the orginal sequence. 
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Figure 3.19--Superimposed power spectral density plots of the original 256-point test sequence generated by 
(3.30) and the high-frequency and low-frequency channel expander outputs. 


a 


Figure 3.19 includes spectral density plots of the orginal sequence and the output of each 
of the expanders. As indicated by (3.10b), the expander reintroduces the original sampling 


frequency at the expense of creation of aliasing terms. For the low-frequency channel, spectral 
peaks occur at the original location, 77/20 (= 0.175 f,), and at 332/20 which, because of aliasing, 
appears as 137/20 (= 0.325 f,). Similarly, for the high-frequency channel, in addition to the 
peak at its original position of 172/20 (= 0.425 f,), an aliased peak is also present at 32/20 (= 


0.075 f.). Although, at the expander output, the width of each peak has been restored to that of 
the input, each peak is six decibels below the original since the signal energy has been divided 


between the original term and an aliasing term. ' 
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Figure 3.20--Superimposed plots of power spectral densities of original sequence S(e’) and of low-frequency 
QMF channel component and high-frequency QMF channel component . 


Superimposed plots of the power spectral densities of the original 256-point test sequence, and 
of the sequences at the outputs of the low- and high-frequency channel synthesis filters are 
exhibited in Figure 3.20. For each channel, the effects of the synthesis filters are evident. The 


aliasing terms have not been completely eliminated, but the spectral peaks of the aliasing terms 
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have been reduced in magnitude by approximately 95 dB. The restored spectral peaks, still 3 dB 


below the original, track quite closely within the passbands of their respective filters. The 


time-domain of the reconstructed sequence, $(n) =So(n)+$,(n), 1s plotted in Figure 3.21. 
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Figure 3.21--Plot of reconstructed version of 256-point test sequence applied to QMF bank input. 


In this plot, the resemblance of the reconstructed signal to the original 1s evident. The 
reconstructed sequence lags the original sequence by 57, the order of the filters from which the 


QMF bank was constructed. Furthermore, if the normalized, mean-squared error is defined as 


= neo =) ; 
ne \é : 
4 | s(n j—s(n) | by fsin}-sin } 
e= s / id n=O 
s=(n)  . 

> s(n) 


n=0 
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then for the sequences plotted in Figure 3.21, €y3 =—-168.2 dB. 
D. M-CHANNEL FILTER BANKS 


Having, in the previous section, discussed and demonstrated the implementation of a 
two-channel, quadrature mirror filter bank structure, it remains to extend the results to structures 
consisting of an arbitrary number of channels. With some simple modifications to relationships 
from sections B and C, M-channel filter banks can be shown to be largely analogous to QMF 


structures. Finally, it 1s worth observing that in literature [13], M-channel filter banks are 


extensively referred to as M-Channel Quadrature Mirror Filter Banks. The use of the term 
quadrature represents, in this case, a misnomer. Nevertheless, the terminology has continued to 


be applied to these more complex systems. 


The structure for an M-channel filter bank is illustrated in Figure 3.22. The structure is 
entirely analogous to the two-channel case. Each channel contains a factor-of-M decimator and 
an expander which, respectively, subsamples by a factor of M and inserts M-1 zeros between 


each sample. The filters are all M"-band filters with, for the ideal case, frequency responses of 


| axe) 5 loatei &-ksos%-(k+1] | 3.33) 


— otherwise 


Furthermore, the concepts of perfect reconstruction and power symmetry also apply. 
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Figure 3.22--Structure of M-channel filter bank. After [11]. 
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Beginning Z-transform analysis, transformation of a decimated sequence requires the 
construction of a "sifting" function equivalent to that defined in (3.5). For M channels this 


function is defined as [12] 


oe sw, (3.34 
Tain) = M M : ) 
where Wy =e!*™“. Employing 3.34 produces s) y(n) = s(M-n)-)y(M'n), thereby ensuring 
that s,,,(n) is zero for non-integer values of n. For any integer n which is an integer multiple of 
M, W,,=1. Otherwise, the series term in (3.34) becomes a summation around the unit circle in 
the complex plane which is evaluated as zero. Evaluating the Z-transform of the factor-of-M 


decimated sequence s(n) produces 


il 
mM 


Sim (Z) s(M-n)-Nm(M-n)- z™ 


(3.35) 


i 
z<|- > 
<uMr 
5M 
an 
= 
a) 
= 
3 
NX 
2 
Ne” 
3 


Furthermore, the Z-transform of the decimated sequence s,,,(n) after expansion by a 


factor of M 1s evaluated as 


[SumJiyy(2) = Zsum(n/M)z7 
é WE s(we . (3.36) 
7 M n=0 M2 
Evaluating 3.35 on the unit circle produces 
| Mal | | 
Sie a= a x S(ei2emm gion | (33% 
™ m=0 


Consequently, decimation spreads, by a factor of M, the power spectral density of a sequence. 


Furthermore, the M-fold decimator produces, from a sequence with a 2-1-periodic frequency 


4] 


response, a sequence with a 2‘M-n-periodic sequence. Similarly, on the unit circle, (3.36) 


becomes 


Mz ga? 
[SuJru( e | 4 ES(eiecinm | 
Bro feientants | 


l 
M m=) 


(3.38) 


By (3.38), it is apparent that expansion of a decimated sequence produces M-1 aliased terms in 


even Intervals across the spectrum. Additionally, the frequency response of the sequence has 


been restored to its 2:%-periodocity. 


Applying (3.35) and (3.36) to the M-channel structure in Figure 3.22 produces the 


following matrix formulation for the system output X(z) [12]: 


Ho(z) FAC) oe Hy-1 (2) 
Xe)= 4 X(2) X(Wy2) = xOWM 2) ] pene) see 7 Hu (Wae2) 
HoCWM 2) Hwt!2) Hy (wet!) 
F 9(Z) 
F | (Z) 
Fyi(z) 


(3.39) 
The matrix H(z) is (3.39) is the M-channel analogy to the alias cancellation matrix. The column 


vector containing the synthesis filter characteristic polynomials can be represented as 


f(z)=[ Baz) Lae Fata For alias cancellation, it is necessary that 
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Zz) are) tz). (3.40) 


Furthermore, analogous to (3.24), it is possible to define a polyphase representation of 
the characteristic polynomial H,(z): 


H,(z) = ~hy(n) Zan 
= h,(Mn) eS z'2hy(Mn + 1Iyz™Mo4 4+ zh (Mn +M-—1)z-Ma 
(3.41) 


M-1 
= Yz™ oh.(Mn+m)z™ 
m=0 n 


SE 2 E, (2) 
m=0 
where Ex m(Z) © ex m(n)=hy(Mn +m). From the definition of the polyphase representation, 


the polyphase representation follows that 


Ho(z) Bie, eeonie ) 2 Eomai(z™) | 
h(z) ea a E1o(z") pe) as Eimi(Z™) i 
z , (3.42) 
Fix (2) Bere Veewem(Z) --- Eyeyitz™) Zane 
= E(z")e@ 
where 
Eye) Eo oe (z™) 
E(z™)= E10(2™) Ei(2™) - Eimi(Z™) 
Eye eEwen(2!)--- Evexnwetiz 
ande(z)={1 z' - meee | 


Finally, the alias compensation matrix H(z) and the polyphase matrix E(z") are linked by 


a relationship analogous to (3.27). First, 
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H(z) = | nz) h(Waz)  nCOWMZ) | 
: ait (3.43) 
= | E(z)e(z) E( (W wz)™ } eCW wz) E((we 2) ) ewe Z) 
.\M | Mo 
Now, since, for any integer k, (wi = Ce =e!?™k =] | (3.43) becomes: 
H™(z) =E(2™)| e(z) e(Waz) -- e(WM2) |, ( 3.44) 
The block matrix on the nght-hand side of (3.44) is equivalent to 
ee ae Wa Wa oo Wy 
ae we? Ww WN 
[ (2) e(Wu2) -- ewittay | = | ft 7 ee 
0 QO. ztMeD we, wv! ee, woeu 
Combining (3.44) and (3.45), therefore, produces a general relationship between the alias 
cancellation and polyphase matrices: 
H(z) = W, Du(z) E7(z™) (3.46) 
where, 
1 QO =: 0 
a ae 
Du(z)= : “ : 
0 0 we ge 


Obviously, from (3.46), if H(z) is paraunitary, then E(z™) is also paraunitary. Additionally, the 
properties expressed by (3.24), (3.25), and (3.29) are equally valid for the M-channel filter bank 


systems. 
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Figure 3.23--Time plot of 256-point test sequence generated using (3.47). 
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Figure 3.24--Power spectral density of 256-point test sequence generated using (3.47). 


To demonstrate the operation of a three-channel filter bank system conforming to the 


structure of Figure 3.22, a 256-point test sequence similar to that of (3.30) was generated: 
s(n) = 1 a a J @{n-94)21 | g2n75 [ cos( 2% 27-n)++cos(s= $5-n)++cos(= 113. n) |. (3.47) 


Equation (3.47) employs the same envelope as (3.30) applied to the sum of three harmonic 


components. Spectral peaks in (3.47) are located at digital frequencies of 2:%-27/256 (=0.027:f.), 
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2-2-55/256 (=0.215-f,) and 2:2: 113/256 (=0.44-f,). Consequently, spectral peaks occur within 


each third of the frequency spectrum below the Nyquist frequency. A time plot of the test 


sequence is displayed in Figure 3.23 while Figure 3.24 presents the power spectral density of the 


test sequence. 
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Figure 3.25--Superimposed plots of impulse responses of filters used to implement three-channel filter bank. 
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Figure 3.26--Equivalent impulse response for three-channel filter bank structure whose filter impulse responses 


are plotted in Figure 3.25. 
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The filter bank employed for the signal decomposition and reconstruction was 
implemented using 14-point third-band filters whose coefficients are tabulated in [11]. The filter 
impulse responses are plotted in Figure 3.25. Satisfying (3.34) for the M-channel case 1s 
significantly more difficult than for the two-channel case. Consequently, as illustrated in Figure 
3.26, the deviation of the equivalent system impulse response deviates from the ideal case of 
(3.31) by an order of magnitude more than was observed for the two-channel structure 
demonstrated in Section C. For the two-channel structure, the root-mean-square deviation from 
(3.31) was 5.16 x 10°. For the three-channel system whose equivalent response is plotted in 
Figure 3.26, the root-mean-square deviation from (3.31) is 7.42 10°. Furthermore, the peak 
amplitude distortion is also greater. The peak amplitude distortion is approximately 4 x 10° dB 
from the power complementary case. This represents a noteworthy increase over the 1 x 10° 


dB peak distortion for the two-channel structure demonstrated in Section C. 


8.004 
a Hpylez 

& ‘ (Z2!H (ers) 1) 8.882 

wo 
a CO 
= 8 @ 
“w v0 
— a - Ss 

is | -@.082 0 == 
= is 


' 
4 


( 


-8. 804 





-8. 806 
8 8.85 Q.1 8.15 0.2 @.25 3.3 8.3 Q.4 @.45 @.5 


Frequency (multiple of fl 


Figure 3.27--Superimposed plots of frequency responses for filters whose impulse responses are plotted in Figure 
3.25 (right-hand axis) and amplitude distortion from three-channel filter bank constructed from filters of Figure 
(3.25) (left-hand axis). 
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Insight is obtained from considering the decimation operation from both the time and 
frequency domains. Time-domain plots for the output of each of the decimators appear in 
Figures 3.28a, b, and c. The power spectral densities of the decimator outputs are superimposed 
in a plot presented in Figure 3.29. As discussed previously, the decimation operator, in the case 
of a factor-of-three decimation, retains only one sample out of every three. Therefore, the length 
of the decimated sequences is one third of the length of the original sequence. Furthermore, the 
effective sampling frequency f,,, for the sequence decimated by a factor of three 1s related to the 
sampling frequency f. by 

ii 
Additionally, decimation produces aliasing terms. However, these terms are of no consequence 


since the lie outside of the region [0, 2-m]. Figure 3.29 displays the power spectral density of the 
content of each filter bank channel with respect to its post-decimation sampling frequency f,,,. 
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Figure 3.28a--Low-frequency channel decimator output for 256-point test sequence generated by (3.47) applied 
three-channel filter bank of Figure 3.25. 
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Figure 3.28b--Medium-frequency channel decimator output for 256-point test sequence generated by (3.47) 
applied to three-channel filter bank of Figure 3.25. 
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Figure 3.28c--High-frequency channel decimator output for 256-point test sequence generated by (3.47) applied to 
three-channel filter bank of Figure 3.25. 


In the case of the low-frequency channel, the analysis filter transmits only that portion of 
the power spectral density of the test sequence which lies in the region [0, 7/3]. As discussed 
during the development of (3.36) and (3.37), because of decimation, the content of the region 


[O, 2/3] will be linearly redistributed over the region [0, 2]. The spectral peak passed through 


the low-frequency channel is, as a result, translated from its original location at 2-m-27/256 


49 


(= 0.105-f,) to an apparent, post-decimation location of 2:1-81/256 (= 0.316-f,,,). Similarly, the 
spectral peak contained in the medium-frequency channel is translated from its original location 
at 2:%:55/256 (= 0.215-f,) to an apparent, post-decimation location of 2--165/256 (= 0.645-f.,,). 


This expectation is confirmed in Figure 3.29 in which the spectral peak of the test-sequence 


component contained medium-frequency channel appears at the predicted location with an image 


appearing at 0.355-f,),. Finally, the test-sequence component passed through the high-frequency 
channel, originally appearing at a location of 2:%-113/256 (= 0.414-f,), after decimation, assumes 
an apparent position of 2:7-339/256 (= 1.324-f,,,). Because of aliasing, this component, in 


Figure 3.29, is indicated at 2-%:339/256 - 2-%-256/256 = 2:%-83/256 (= 0.324-f,,,). 
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Figure 3.29--Superimposed plots of the power spectral densities of the original sequence and the outputs of each 
decimator for the 256-point test sequence of (3.47) applied to the three-channel filter bank of Figure 3.25. 


In Figures 3.28a, b and c, the zero-crossings of the sequences appear to occur at similar 
frequencies. This observation is confirmed by plots of the decimator output spectral densities 


superimposed in Figure 3.29. In fact, for the case under consideration, the spectral peaks of the 
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decimator outputs are separated by a maximum of 0.5-f,. This example provides insight into the 
nature of the decimation operator. In general, decimation transforms a narrowband process such 


that the result "fills" the spectrum below the Nyquist frequency of the channel. 
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Figure 3.30--Superimposed plots of the power spectral densities of the original sequence and the expander outputs 
for each channel of the three-channel filter bank of Figure 3.25. 


The power spectral densities for the expander outputs for each channel of the filter bank 
represented by Figure 3.25 are plotted in Figure 3.30 above. As predicted by (3.38), application 
of the expansion operator to a decimated sequence imposes two consequences. First, the 
effective sampling frequency of a decimated sequence is related to the sampling frequency of the 
original sequence by 

[fausdts = 3 faus = f 
in the case of expansion by a factor of three. Secondly, the aliasing terms generated by 


decimation, which are originally outside of the region [0, 2-7], are translated to within the region 


Ores | 
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In the case of the test sequence generated by (3.74), two aliasing components in addition 


to the desired component have appeared in each channel. For each channel, the aliasing 
components appear at integer translates of 2:7/3 with respect to the restored component of the 
original sequence. Because the spectral peaks of the original signal were separated from each 
other by roughly 7/3, each aliasing term generated for each component coincides fairly closely 


with one of the other components. This occurrence 1s reflected in Figure 3.30. 

In the low-frequency channel, the location of the peak of the spectral component has 
been restored to its original location of 2:%-27/256=2:%-81/768. However it 1s accompanied by 
aliasing terms at 2:%-337/768 (=0.438-f.) and 2-1-593/768 (=0.772-f,.) whose image appears at 
2:-m — 2:7°593/768 = 2:-m-175/768. The restored medium-frequency-channel component, which 
reappears at 2:%-55/256=2-m- 165/768, is accompanied by aliasing terms at 2-:%7-421/768 
(=0.548-f.) whose image 1s present at 2:%-347/768 (=0.452-f,) and at 2:%-677/768 (=0.882:f,) 
which has an image at 2°7-91/768 (=0.1185-f,). Finally, to the restored component for the 
high-frequency channel which is located at 2:m-113/256=2-m-339/768, are added aliasing terms at 
2:%°595/768 (=0.775:f.) for which an image appears at 2:%-173/768 (=0.225-f,) and at 2-7-83/768 
(=0.1098-f,). 

Figures 3.31a superimposes plots of the outputs of the synthesis filters for each channel 


of the filter bank of Figure 3.25. Each plot indicates the spectral content of the corresponding 


channel. The results of the recombination of the channels are plotted in Figure 3.31b. In Figure 


D2 


3.31a, it is apparent that the aliasing components of the expander outputs have not been 
completely blocked by the synthesis filters. In the worst case, for the spectral region containing 
the medium-frequency channel, residual energy from an aliasing component is only 
approximately 20 dB below the desired spectral peak for that channel. However, upon 
examination of Figure 3.31b, it becomes apparent that alias cancellation does occur. The power 


spectral density of the reconstructed sequence very closely coincides with that of the original 
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Figure 3.31a--Superimposed plots of the power spectral densities of the original sequence and of the expander 
outputs for each channel of the three-channel filter bank of Figure 3.25. 


Figure 3.32 presents a time plot of the reconstructed sequence. Again, the reconstructed 
sequence appears to be an approximate delay of the original sequence. In fact, when the signals 
are synchronized, the normalized mean-square error (3.32) of the reconstructed signal is -66.94 
dB. The the reconstruction error exceeds that of the two-channel demonstration of Section C 
because of the poorer quality of the filters with which the three-channel filter bank of Section D 


has been implemented. As indicated by Figure 3.26, the equivalent impulse response of the 


>3 


three-channel filter bank deviates from a pure delay by a margin three orders of magnitude 


greater than the deviation for the two-channel filter bank. 
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Figure 3.31b--Superimposed plots of power spectral densities of original sequence generated by (3.47) and of 
Sequence reconstructed by thee-channel filter bank of Figure 3.25. 
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IV. THE THEORY OF MULTIRESOLUTION SIGNAL PROCESSING 
A. INTRODUCTION 


Having, in the preceding two chapters, laid the necessary groundwork, the theory of 
multiresolution signal processing will next be considered. In this chapter, Section B, presents 
the concept of multiresolution analysis. In Section C, Mallat's multiresolution algorithm will be 
developed from a projection operator perspective and the equivalence of multiresolution 
mathematical operations and two-channel QMF banks will be demonstrated. Section D will 
outline the development of the Laplacian pyramid and the A Trous algorithm, two of the earliest 
multiresolution decomposition techniques. In Seeiion E, multiresolution structures comprised of 


cascades of filter banks will be constructed and demonstrated. 


Signal processing techniques commonly entail decomposing a signal by representing it in 
terms of its projection on a vector space. The most common method, the Fourier Transform, 
defined in Chapter II: 

fw) = A {£t)} =J7 F(t) eI?" de. (4.1) 
In analyzing a time-varying signal, however, (4.1) presents an obvious disadvantage: Onlv one 
representation vector is used for all time. Consequently, time-varying aspects of f(t) are 
averaged over all time and lost. To address this shortcoming, the concept of the Shovt-Tinte 
Fourier Transform (STFT) was developed [14]: 


F(o, D=f7 f(t) w(t—t)e7I?! dt. (4.2) 


tn 
ta 


Equation (4.2), through the introduction of a window function w(t), improves on (4.1). 
Typically, the window a function is either strictly time-limited or possesses a rate of decay such 
that its value outside of a limited, contiguous region is negligible. Additionally, equation (4.2) 
can be interpreted as the projection of the signal of interest on modulations and translations of a 
vector w(t). Employing concepts from Chapter II, (4.2) can be interpreted as a mapping from a 


one-dimensional space of real numbers to a two-dimensional space of real numbers, or, 
F: R > R’ 
Additionally, the representation vector r(@, Tt) for the projection operator F, the operation of 
(4.2) 1s 
r(@, t)-= w(t- tT) 3?" (4.3) 
In the branch of mathematics known as group theory, the operation (4.2) belongs to a 
particular class of operators known 1s the Weyl-Heisenberg Group [15]. A group is a set of 
transformations satisfying the properties of closure, associativity, identity and invertibility [16, 
17]. The Weyl-Heisenberg Group consists of a family of transform operators characterized by 


modulation and translation of a single representation vector [18]. 


Although it is commonly employed, two shortcomings of the STFT have been asserted. 
The representation vectors for (4.2) represent a frame in the sense described in Chapter II [6]. 
Therefore, if it is necessary invert the transform in order to reconstruct a signal from its STFT 
decomposition, a dual operator must be constructed. Secondly, a sampled, discretized STFT 
Operator partitions an analyzed function's two-dimensional conjugate space into uniform, 


rectangular partitions [19]. In the conjugate space, the spectral bin partition dimensions are 
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inflexibly dependent upon the window function w(t) and to not vary with either translation or 
modulation. Many processes (for example, biological processes) can be characterized by 
components whose bandwidths increase with frequency. Consequently, (4.2) provides 
representation which 1s less than optimum processes comprised of spectral components of 


varying bandwidths [20]. 


For the reasons described above, alternative methods of signal decomposition have been 
suggested. Instead of a representation vector of the form (4.3), employment of a representation 


based on a family of functions [21] 





vao(t) = fal"? y( 432 | (4.4) 
produces a transformation W: R > R’ such that 
WhE)Y = Wa, by='lal-#2(7_ f(t) v (2 ] dt. es) 


Similar to the STFT, transformations of the form of (4.5) also comprise a distinct class in the 
field of group theory. Transformations based on scaling and translation of a common 


representation vector comprise the affine group [18]. 


If, as in the case of digital signal processing, it 1s desired to restrict the transform to a 
lattice of discrete points, the representation vector becomes 
Winalt) =a” wap t—n bo). (4.6) 
The representation (4.6) results in a transformation W: R > Z’, where Z is the set of all 
integers, such that 


W{f(t)} = wm, n) =a5™* J@ f(t) y*(ag™ t—n bo) dt. (4.7) 
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Additionally, (4.6) constitutes a version of (4.4) sampled at a = a," for a, > 1 and at b = 


n'by'a)” for b) #0 . To conform to the conventions of octave-band filtering, a, is typically 


selected as a,=2. The selection of b, determines whether {Von atne z constitutes an incomplete 


set, a complete orthogonal set, or a frame. The representations of greatest interest to 
multiresolution signal processing are chosen such that b, = 1, providing unit translations with 


respect to sampled data. 


The time-frequency properties of the representation vector y_, , address some of the 
shortcomings sometimes ascribed to the STFT representation. As the scaling integer m 
increases, the representation y_ , becomes more and more spread in the time domain, and, 


consequently, more concentrated in the frequency domain. Consequently, projection on a highly 
dilated vector function will provide poor time resolution but sharp frequency resolution. 
Decreasing m causes the reverse effect: Concentration in time and spreading in the frequency 
domain. In the case of a highly contracted representation vector, the transform operator provides 


sharp temporal resolution but poorer spectral resolution. Furthermore, the spreading effect of 


the representation y,, , occurs in a logarithmic manner. The bandwidth of the representation 


Winn Will be proportional to its center frequency. 
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Vin 


Figure ¢./--Venn diagram illustration of concept of embedded vector spaces. After [24]. 


B. THEORY OF MULTIRESOLUTION ANALYSIS 


In order to characterize the vector spaces consisting of the span of {y,, ,$., .-7: it is first 
necessary to consider another set of basis functions ia which spans another set of vector 


spaces {V_}.., € L?(R). The operator A,_, is defined as the projection of some function f(t) on 


V : 


m* 


a 
Am-1{ f(t) } ale dma )Omalt) (4.8) 
To be a multiresolution analysis, the set of operators {A,} must satisfy six properties [22,25]: 


l. A_, is a linear operator which uniquely and completely approximates f(t) at a resolution 


in-1 


of 2™. Consequently, the approximation 
eee) mene) (4.9) 
In words, A, ,{f(t)} contains all of the information about f(t) which can be obtained at a 


resolution 2". Repeated projection upon V,, does not add or subtract any information to 


Awl (f(t) } : 


bo 


Normalized Magni tude 


Of all possible functions which exist at a resolution of 2”, A__,{f(t)} is the function 


which most closely resembles f(t): 


V g(t)e Vax |[g(t) - f(t)|] = JA, {£0t)} - £(¢)]] (4.10) 


Approximation of some signal f(t) at one resolution 2™ contains all information 
necessary to approximate it at the next resolution 2™”'. This concept suggests a family of 
embedded, closed subspaces: 


V..,cV éL(R) V me Z, (4.11) 


m~+t 
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Figure 4.2--Spectral illustration of concept of embedded vector spaces. The vector spaces represented are 
based on Daubechies' orthonormal scaling function on (0, 11]. After [25]. 


Illustration of the concept of embedded spaces can be accomplished by either of two 


methods. The first illustration is via Venn diagram. In Figure 4.1, the embedded ellipses 


illustrate two related vector spaces. The outer ellipse represents the span of the vector space V,, 


and the inner ellipse the span of vector space V_,,. As indicated by the diagram, the information 


which can be extracted from projection on vector space V,,,, is less than what can be extracted 
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from projection on vector space V_,. The two spaces differ by the information lost when some 
function f(t) is approximated at 2™”' instead of at 2". The vector space V_ possesses a greater 
density than the vector space V_,,,. Consequently, the approximation at 2 will contain a greater 


quantity of information about the original function than the approximation at 2™”', 


Figure 4.2 provides an illustration of the concept of embedded subspaces from the 
perspective of partitioning of the frequency spectrum. The lower half-band represents the vector 
space V_, while the lower fourth-band represents V___,. Approximation A,__,{f(t)} at resolution 
2™, therefore, entails an approximation based on the spectral components of f(t) contained in the 
lower half of the frequency spectrum below the Nyquist frequency. Similarly, approximation at 
a resolution 2™”' entails a representation based on the lower fourth of the frequency spectrum 
below the Nyquist frequency. Consequently, an approximation of f(t) based on Figure 4.2 at 


resolution 2™ contains only the spectral content of f(t) in the range [0, m/2]. An approximation 


m-+] 


of f(t) at resolution 2"” contains the only the spectral content of f(t) in the range [0, 7/4]. 


4, The approximation operation 1s similar at all resolutions. The spaces of approximated 
functions can, therefore, be derived from one another by scaling each approximated 
function: 

Vmed, f(theV, = f(2-t)e V,.. (4.12) 
5s The approximation A,_,{f(t)} of a signal can be characterized by 2™ samples per unit 


interval. When f(t) is translated by an amount proportional to 2™, A__,{f(t)} 1s translated 
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by the same amount and Is characterized by translations of the vector space projections. 
More simply, 
Vne Zz f(the V, & f(t-2™n)e V,. (4.13) 
Equations (4.12) and (4.13) suggest a family of basis functions ,_(t) similar to that 
characterized by (4.6): 
Oma(t)=2-™2O(2™t-n) V m,neZ’. (4.14) 
The vector space oaeae z€ V,, consists of integer translates of >. (The t-dependence has 


been suppressed for compactness of notation.) The approximation operator A__,{f(t)} is simply 
a projection of f(t) in the space of the vectors {@, ,}. Additionally, (4.12) reinforces the concept 
illustrated in Figure 4.2. Equation (3.7) indicates that a time-domain contraction of a signal 
causes a dilation--or spreading--of that signal's frequency spectrum. The frequency spectrum of 
V 4.» therefore, occupies half of the bandwidth of the frequency spectrum of V_. Furthermore, 
if the space spanned by V_ coincides with a lowpass region of the frequency spectrum, V_, 


will also coincide with a lowpass region. 


Finally, one additional property remains to complete the definition of a multiresolution 


analysis: 


6. A continuous function f(t) can be initially considered to be represented with infinite 
resolution. Regardless of the scale, all information regarding f(t) 1s onginally assumed to 


be known. Applying the approximation A, {f(t)} results in some loss of information. 
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Increasing the degree of coarseness of A, increases the amount of lost information. 


Consequently, 


lim Vm= U_ Ve =L*R) (4.15a) 
and 
lim Vm= A Vi={ 0 f. (4.15b) 


Equations (4.15a) and (4.15b) follow, by induction, directly as consequences of (4.11). 
Equation (4.15a) states that, because of (4.11), the manner in which ea 7 are related to 


each other, lim Vm consists of the union of the spans of all {v_} By induction from 


Me 2: 


(4.11), all vector spaces {V_}_ <z are subsets of lim Vm. Using the notation of set theory, 
the concept behind (4.15a) 1s expressed as 

\Vingmez © ,lim Vm. 
Consequently, adding to (4.15a) the concept expressed by (4.12), as the resolution of an 


approximation A__,, the projection on vector space V_, becomes infinite, it is possible to 


m-1° 
represent f(t) with arbitrary precision. Furthermore, employing a concept from Chapter II, 


it is possible to find another 


mp? 


lim_ Vm is dense in L*(R). In other words, given any vector 


vector o which is arbitrarily close to >, ,. Contained within the union of all definable vector 


m,kz#no 


spaces {V_} is all the information which is known about a function. 


The indication of equation (4.15b) exactly opposite of that of (4.15a). By induction 


from (4.11), tim Vm constitutes the /east common subset of all of the vector spaces 


{Vin} me z-Consequently, as the coarseness of the approximation A, ,(t)} becomes infinite, all 
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information about f(t) is lost. From the perspective of Figure 4.2, the portion of the frequency 


co 2m 


spectrum spanned by tim Vm is contained in the closed interval | 0, lin oa | 


Consequently, at the limit, the approximation of f(t) is characterized by only its DC component, 


or equivalently, lim Am{f(t)} becomes a constant-valued function. 


Regarding multiresolution analyses, Mallat [22] proved a theorem which provides the 


theoretical foundation for all further development. The theorem states that, given a 


multiresolution approximation of L*(R) by projection on vector space V,,, there exists a unique 
function o(x)€ L*(R), called a scaling function, such that 1Onn(X) }ne7 as defined by (4.14) 
constitutes an orthonormal basis of V_. That ts, 


C= Omk) =Sks V n, kom (4.16) 


Furthermore, V_ consists of the closed, linear span of Gn tas 


=. —o {t) 
ote” 4h -— 0% 4 (t) 
@ 


' =I 


Scaling function value 


2 4 6 8 12 
Scaling function domain, t 


Figure 4.3--Superimposed plots of Daubechies' orthonomnal scaling function 9, ,(t) supported on [0, 11] and 
of ©, o(t). 
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To characterize the relationship between the scaling function spaces {V_}, it is useful to 
examine the projections of their vectors onto each other. From (4.11), (4.12) and (4.13), it is 
evident that 

Omernalt)€ Vin SS On Athe V,, (4.17) 


or equivalently, les within the span of {@, ,}. Figure 4.3 provides an illustration of a 


m-l.na 
scaling function at two adjacent scales. The example presented is based on Daubechies' 
orthonormal scaling function supported on [0, 11]. 9., o(t) 1s obviously a contracted version of 


 o(t) with its amplitude increased by a factor of nD . Because of (4.17), the Fourier series 


theorem (2.17) can be applied to obtain 


Om+1.n{t) =E(Onsin mk J Oma) (4.18) 
k 


Substituting (4.14), the definition of O__, into the inner product term of the summation (4.18) 


m,n? 


and applying a change of variables to evaluate the resulting integral produces [22] 
(Gmei.a: Omi ==: (1.0, boacza } (4.19) 


Equation (4.19) clearly indicates that the coefficients of the series (4.18) are independent of 


scale m. The summation (4.18), therefore, becomes 
Om+i.n(t) = F E( 1.0, barra | Om, k(t). (4.20) 
k 
Next, substituting (4.14) for the appropriate terms on each side of (4.20) results in 


ESO) (Rie x (O10, borcza | -o(2-™ t—k). (4.21) 


Applying the change of variables u=2-". t-2-n and the translation of indices k=k'-2-n to (4.21) 
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finally yields 
o(u) =E (10 do.i-20 O(2U-K) (4.22) 
For later convenience, (4.22) can be rewritten as 


> O($) = = h(k) O(u-k) where hWEH{ O10, Do. k ) (4.23) 


Selecting this normalization for h(k) forces ~ h(k) =1. Asa result, 


| Eh({k) ei 





w=0 . : H(e!®) ae = 


Furthermore, since (u) is real, {h(k)}, . 7 is real. And finally, if o(u) is compactly supported, 


then the filter h(k) is a FIR filter. This occurs since, for compactly supported 9, in the inner 


product term of (4.23), there will be only a finite number of translations which will be evaluated 


as non-zero. The details of the development of (4.19) - (4.23) are presented in Appendix A. 


Equations of the form (4.23) comprise in literature a class of equations referred to as 
two-scale difference equations [26] or as dilation equations [27]. This class of equations will be 
considered in greater detail in Chapter V within the context of basis functions for multiresolution 
analyses. From the preceding development it may be concluded that the functions @ which 
satisfy the properties for a multiresolution analysis consist of solutions to two-scale difference 


equations. 


To obtain a multiresolution transform it is necessary to combine (4.20) with (4.23) to 


obtain 


66 


Ometn(t) = V2 Lh(kK- 20) Ome(t). (4.24) 
Applying Parseval's equality (2.18) to (4.24) results in a relationship between the projections of 
some function f(t) onto V_ and V_,, [22]: 


, met a si) ~ h(k 215)) G Omk). (4.25) 
By (4.25), given A,_,{f}, the determination of the coefficient (f, 0,, ,) for each term of the 


series for A, {f}, requires filtering and decimation of the sequence {(f, On Oh: The block 


diagram for the operation (4.25) is depicted in Figure 4.4. The similarity between Figure 4.4 


and a channel of a QMF analysis bank 1s evident. 


{( 6 Om neo hin) -e|2t p> (£0... JF > 


Figure 4.4--Block diagram of multiresolution transform characterized by (4.25). 


Now, to complete the development of the concept of multiresolution analysis it is first 
necessary to consider some of the consequences of the preceding development. First, consider 
the concept of embedded subspaces characterized by (4.11). Let W_ be the orthogonal 


complement of V_ such that 


V,AW_=2 (4.26a) 
and 
Vagere— Vv (4+.26b) 


From (2.26b), by induction, any scaling function vector space V, can be decomposed as [23] 


V7 = ) Ay in: (4.27) 
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And finally, the direct sum of all definable {W,,},,.7 is dense in L’(R) [21]: 
@. Wm =L(R). (4.28) 
Figure 4.5 presents a two-scale example of (4.26b) from the perspective of spectral 

contents of V_ and W,,. The complementary subspace W,, is defined such that it is closed and is 
spanned by a set of vectors {W,, ,}ne7. Members of the set of vectors {y,, ,} are related to each 
other in a manner similar to (4.14): 

Wina(t)=2-™2 y(2-™t-n) V m,neZ?’ (4.29) 
Overlap exists between the spectral regions characterizing V_., and W_.,. This overlap occurs 
since the plots consist of frequency plots of FIR filters. However the FIR filters represented are 


related to orthogonal basis function and consequently, the filters for V_., and W___, are 


themselves orthogonal to each other in the time domain. 


aoe 
ws? 


—Spectral cantent of vector space V 
—Spectral content of vector space V+1" 
Spectral cantent af vectar space W+l 







8.8 





8.6 


a.4 


Normalized Magnitude 


a.2 


Q 8.95 Q.1 8.15 Q.2 8.25 Q.3 a.35 a4 Q.45 a.5 
Digital frequency, (multiple of fl 


Figure 4.5--Spectral content of example of concept of orthogonally complementary, embedded subspaces. The 
example is based upon Daubechies' wavelet and scaling function on (0, 11]. The lowpass processes are 
>, , and o,,.,, and the bandpass process Is Y,,.,., 
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The set of vectors {W,, }nez Constitutes the set of wavelet functions for the 
multiresolution analysis. From the relationships between the spaces W,, of wavelet vectors y,, . 


and the spaces V__ of scaling function vectors >, , a number of relationships between the sets of 


mn 
vectors o>, , and y,, , may be deduced. Between the scaling and wavelet function vectors, 
follows from (4.26a), the relationship exists: 

(G:.0 Wm.) = 0. (4.30) 
Within a common scale, therefore, wavelets and scaling functions are orthogonal to each other. 


Furthermore, as a consequence of (4.11) and (4.26b), 
(Vas Vie) = 8,m Be (4.31) 
In words, unlike scaling functions which are only orthogonal to each other with respect to 
translation within a common vector space V_, wavelets are orthogonal to each other with respect 
to translation and scale. Wavelets not contained within the same vector space W,, are by 
definition orthogonal to each other. Finally, because y, .€ W, CV. 
(0, Via) #0, V ksm-l. 

More generally, if the wavelet y, ; is supported on the interval (0, L] e R, then, 

(Ono Ve, #0, V kml, In-jl <L. (4,32) 


Wavelet vectors lie within the span of higher resolution spaces of scaling functions. 
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The projection of some signal f(t) on the wavelet space W_ constitutes the detail signal 
[22] and contains the difference of the information contained in A, ,{f} and that in A_{f}. The 


detail signal at resolution m is denoted by 
Dri {fl)}25 (f, Ya ) Walt) (4.33) 
Equation (4.26b) implies that [21] 
Daa) (4.34a) 
or equivalently, 


2 (Want.a | WV mela = 2 Ca | Omen a 2 (Ona. | Orel . (4.34b) 


Substituting (4.24) into (4.34b) and expanding the result produces a time-domain form for the 


wavelet decomposition: 
5 (Ves £)Vnetn =F] (Omar f] ma—2EAk-2n) hj-2m) (Ons, F)Oms |. (435) 
n ii 


Obviously, except for simple cases such as the Haar basis (demonstrated in [(21]), (4.35) is not 
easily separable. In order to obtain a relationship between scaling and wavelet functions it 


therefore becomes necessary to examine their Fourier-domain properties. 


Equation (4.22) represents the starting point for considering the frequency-domain 
properties of the scaling function . If the Fourier transform of 6 is defined as in (4.1) 
(cw) = f= o(t)eF* dt, 
evaluating the Fourier transform of each side of (4.23) produces [23] 


(20) = H(ei°) o(@) (4.36b) 
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where the Fourier series H(e’) of h(k) is defined as 
H(e! “)=5 h(k) Jee, (4.37) 
Next, Poisson's summation can be employed to present an alternative expression for the 
orthogonality of ® with respect to translation [23]. This development begins with the complex 
Fourier series expansion for a train of Dirac impulse functions 6(a@)) [28]: 
E5(w-2nk) = 57 Deke, (4.38) 
Poisson's summation involves convolving a function with each side of (4.38). Convolving the 


A 2 ; P 
expression | o(@) | with the right-hand side of (4.38) produces 


"#E5(@-2ek)=E| §@-2n16 | (4.39) 





| o(@) 
, 2 
Furthermore, convolving | () | with the nght-hand side of (4.38) results in the expression 


| Ga) | * Eee = ZS] Fe | eto ak 


, | oe : (4.40) 
X ek fff O(a) o(v) eS du dv e¥5 dé 


Applying the change of variables of integration w = v - u and regrouping the terms contained in 


the second line of (4.40) produces 





| 9(@) * Deis =ZelX® ff d(u) o(u+ w) f eS dE dudw. (4.41) 


Now, since 
f JOOS dE =2n6(w-k), 


by the integral sifting property of the Dirac delta function 5(w), (4.41) becomes 


| Ga) | * Zeke = 2eE ck ff g(uo(u+w) du Sw—K) dw 


2nd ek® f ou) d(ut+k) du 


(4.42) 


71 


Finally, by (4.16), the integral in the final line of (4.42) is simply (0, 9, % ,) = 0), where Cmmmg 
Kronecker's delta function. Consequently, the series in the final line of (4.42) contains only one 


non-zero term: when the index k= 0. Therefore, 
, 2 | 
| o() | «Del = 2. (4,43) 
Substitution of (4.39) and (4.43) back into (4.38) produces 
x 2 
BH o(w- 27k) | =1. (4.44) 
Applying the results of (4.44) to (4.36) yields 
, 2 | oT 2 ; 
x | 6(2@-2nk) | =z | Hel") || O(@ - ak) ae (4.45) 
The right-hand side of (4.45) can be expressed as 


x | Helter") * | o(w - 2kn) +E Ei(ey ot!) al o(@ —(2k+ 1)n) =e 


(4.46) 


Since, as a consequence of its definition, H(e’) is periodic with respect to 2-2, the terms of 





| Hie Orr | “and | H(e® 2k) |? Can be factored out of the series on the left-hand side of 
(4.46) producing 
| 2 ‘ 2 | 2 , 2 
| Hei9) |Z] d(@-2km) | +] He) || d@-2k+De) | =1. G47) 
The remaining two sums in (4.47) are identical to (4.44) and are evaluated as unity. Therefore, 


(4.47) becomes 


| H(ei®) +] H(ei (7) oe (4.48) 


ii: 


Comparing (4.48) with (3.20b) reveals that the frequency response H(e’®) for the 


approximation filter defined by (4.23) 1s power-complementary with a version of itself shifted in 
the frequency domain by m. Therefore, H(e’”) either is constant or possesses a half-band 


frequency response. Furthermore, if the scaling function © is normalized such that 


2 


= |, then, as a consequence of (4.36b), 





| (0) 


| Eice®) | =a, eo) 


a a ER SR ear: SESS || EESRERP S a 


It remains to characterize the wavelet functions. This can be accomplished by 
completing the development of the Be ificesolution decomposition begun in the previous section. 
The first true multiresolution algorithm was developed by Mallat [22]. Two ney igus | 
will be considered tn the next section which were not directly derived from the mathematical 


definition of a multiresolution analysis (4.8)-(4.13). 

The multiresolution decomposition process begins with the projection of a function f(t) 
on the scaling function space of finest resolution {0, Byer [21]. The functional expansion 
appears as 

f(t) = A_,{f(t)} = 2 Con O(t-—n) where conf Ov.n i (4.50) 


Generally, the decomposition is performed on a sampled data sequence and the scale of the 


initial approximation ts defined such that the first equality in (4.50) holds. Since the basis 


ig 


vectors 0, , € Vo, the original approximation (4.50) is also clearly contained within the span of 


ws 


8 
The next step in the multiresolution decomposition requires the division of the scaling function 
space V, into orthogonally complementary subspaces V, ® W, = V,. Re-arranging (4.34a) leads 
to 

A {f(t)} = Agi f(t)} + Dol f(t) }. 
Consistent with the definition (4.8) of the approximation transform A_ and the definition (4.33) 
for the detail transform D,, 


Ao{f} = 21k O1.k 


and 


Do{f} = 2 Dik Wik 


WheDecye = (f, >, J and b,,=(f, WV, ): A relationship between 0, , and 0, , was presented in 


(4.23). Furthermore, (4.32) showed that y, , lies within the span of Vo. Consequently, 
following development analogous to that for (4.23), the Fourier series expression for wavelet 


function becomes, 


7 (5) = Ze(k) O(u-k) where g(k)=4 (Wi.0 ot) (4.5 1a) 
since 
(Yio bm | = = (yi, borcza J (4.51b) 


Since, if y,, , 1s real, the filter coefficients {e(k)}, . , are strictly real. Additionally, Fourier 
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transformation of (4.51a) produces a relationship analogous to (4.36): 
W(2 w) = Gle!®) 0(@). (4.5 1c) 
Finally, the form of (4.24) for the wavelet function appears as 


uewnwna)y = 2 E glk -21) Omx(t). (4,52) 


Applying Parseval’s equality (2.18) to (4.24) and (4.52) in order to calculate c, , and b, , in terms 
of c, , produces 
3) ai Zh(k - 2 N) Cok 


and 


J2 :Zg(k-2n) co, | 


bia 


The approximation based on the expansion of the terms {c, ,°0, ,} constitutes a "blurred" version 
of the expansion of the terms {c, ,"d, ,}. The detail which is lost in the approximation A) {f} is 
contained in an expansion of the terms {b, ,y, ,}. This process can be repeated as many times 
as desired. The resultant general expression becomes 


Cie (2 -Zh(k-2n) cm. 
and (4.53) 
Gaeta (2 -Zg(k-2n) em 


where each successive iteration generates in an approximation of f(t) containing half the 
resolution of the previous approximation. At each iteration, the approximation and detail of the 


signal are related to the previous approximation by 


esi ty eee ei (4.54) 


) 


Furthermore, the operations described by (4.53) closely resemble the analysis bank for the QMF 
bank structure investigated in the previous chapter. The approximation and detail operations at 
each resolution consist of filtering a common sequence with two different filters. The output of 
each filter is then decimated. The resemblance of QMF structure to the operation characterized 
by (4.53) cannot yet be asserted to be exact since the appropriate relationship between the filters 


whose impulse responses are h(k) and g(k) yet to be established. 


To reconstruct a signal from its decomposition, the definition (4.8) of the approximation 
operator is employed. The expansion coefficient c, , is, by definition, 


Cok = (Ana, Ona) = CARS me) + (De @,, ¢): (4.55) 


Substituting into (4.55) the expressions for the series A, and D., produces 
Cm = 2 Cmel.k Co Derik mi 2 Dinetk Coe Vani (4.56) 


The inner product terms in (4.56) are precisely those encountered in (4.19) and in (4.51b) except 
that the positions of the translation indices n and k have been reversed. Consequently, (4.56) is 
equivalent to 


Cma == Lh(n-2k) Cmetk t+ Dg(n-2k) bere (4.57) 
42 k 7 k 


The operation in (4.57) entails expansion of the detail and approximation series coefficients 
followed by the application of filters whose impulse responses are the time reversals of those 
used in (4.53) for signal decomposition. Again, a resemblance to the QMF bank synthesis 


structure 1s noteworthy. 


Examining the decomposition operation in the frequency domain provides further insight 
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into the nature of the wavelet functions [23]. First the Fourier transform Aj(q@) of the 


approximation 


A_(th=A_,{f(t)}=f(t) is expressed by 


A_,(@) 2 Co.n i Oo.n(t) ceo dt 
Fico, eI" J (thet de (4.58) 


Co(ei®) O(a) 


where C,,(e! JS Cmke?°*. Applying similar analysis, Ao(@) and D(a), the Fourier 
k 


transforms of the approximation A,(t)=A,{f(t)} and D,(t)=D,{f(t)}, respectively, appear as 


Ao(@) 2 C1. Jo. Dik(t) eI! de 

= Eek fr OCS —k) ed?! de 

= J2 Yc,,e72ek JP o(u) e722" du 
k 

J2 Cy (ei?) 9(2 @) 


and 


(4.59a) 


Ao(@) 


Do(w) = J2 Bilei2*) H(2@) - (4.59b) 
where Bm(e!®)=E mi eJ°K By (4.54) therefore, it follows that, 
A-1(@) = Ao(@) + Do(@), 
or equivalently 
f(@) = Co(e!®) 9(@) = (2 Cie?) (2 @) + /2 Bi (e!?*) H(20). (4.60) 
Substitution of (4.36) and (4.51c) into (4.60) produces 


Co(ei®) 6(w) = (2 Ci(ei2®) H(ei*) o(w) + (2 Bi(ei2”) G(e!®) o(@), 


or, after dropping the common factor of 0(), 


Co(ei®) = 2 Ci(ei2®) H(ei®) + 2 By (ei?) Glel®). (4.61a) 


oa 


From their definitions, C,(e®) € L*([0, 2-]) and is periodic with respect to 2-m while 
C,(e *°), B,(e’*°) © L?({0, m]) and are periodic with respect to m. Therefore, it follows from 


(4.6la) that 


; Co(e!®) || H(ei®) Ge) |. aan 


Co (ei (etn) H(e! (@+")) G(e! (car) B,(e!2®) 
From (4.61b), define 


H(ei®) Gel®) | 


a3 
H(e! 4 H(e! (@+")) Gieule)) 


Now, if the H(e’®) is unitary in the sense described in Chapter III, then multiplication of each 
side of (4.61) with its Hermitian transpose yields 
| 2 , 2 , 2 
| Co(ei®) | +] Co(ei@*”) | aN Ci(ei2®) | +] By(ei2®) *]. (4.62) 

Integrating (4.62) over the interval [0, x] produces 

' | 2 | 2 ‘ | 2 : 2 

Jo | Cole?) | +] Cole) [odo =2 Jp | Cr(ei2e) | +| Bicei2®) | do 
which is equivalent to 


Bi(e2%) | dw. (4.63) 





- | 2 | 2 
Jo" | Colei2) | do=2 J> | Ci(ei2®) |° do+2 ff 

Equation (4.63) 1s an expression of the orthogonality of the decomposition: The energy in the 

approximation A,{f(t)} is simply equal to the sum of the energies in the approximation A, { f(t)} 


and the detail D,{f(t)}. No energy is present in terms comprised of the cross-product 
| C (e?*) | | Be”) | . Furthermore, the specification of H(e’®) as unitary is exactly 


equivalent to the specification that the alias compensation matrix for the two-channel QMF bank 
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is exactly equivalent to the specification that the alias compensation matrix for the two-channel 
QMF bank be paraunitary. Consequently, a single-stage decomposition followed by a 


single-stage reconstruction employing Mallat's algorithm 1s exactly equivalent to a QMF bank. 


\foreover, the paraunitary nature of H(e'”) implies 


| Her) | +] Geis) | =1 


| | | (4.64) 
ACetSnG(e 3°) + Hle#?™) Glew) = 0 
One solution for G(e’®) as constrained by (4.64) 1s given by 
Gel®) =e °F H(hO7”) (4.65) 


where L represents the length of the FIR filter h(n). Selection of (4.65) dictates that the filter 


whose impulse response ts g(n) will be a FIR filter of length L and that, by (4.5lc), the wavelet 


function y(t) will have exactly the region of support of the scaling function 0(t). 


D. EARLY MULTIRESOLUTION ALGORITHMS 


{9 





Figure 4.6--Block diagram of process by which detail is extracted from signal through the use of the Laplacian 
pyrainid algoritlun. After [29]. 


In this section, two early multiresolution techniques--the Laplacian pyramid and the a 
trous algorithm--will be considered. The first, the Laplacian pyramid, represents an algorithm 
developed in the early 1980's for the purposes of image coding. Conceptually, the Laplican 


pyramid very closely resembles Mallat's multiresolution algorithm presented in the previous 


is 


section. In fact, Daubechies credits the Laplacian pyramid with providing Mallat with 


significant for the development of his own multiresolution algorithm [21]. 


The Laplacian pyramid, as with Mallat's multiresolution algorithm, entails iterative 
filtering of a sequence to progressively smooth out the rapidly varying components [29]. At any 


stage, the next coarser approximation is obtained from 


fk+i(n) = 2 w(2n ~ m) fx (m). (4.66) 
Unlike Mallat's algorithm, however, where the detail 1s extracted through filtering, the Laplacian 
pyramid extracts the detail through re-expanding and filtering the approximation and subtracting 


the result from the previous approximation. This process 1s illustrated in Figure 4.6. 


Mathematically, extraction of the detail lost A,,,(n) when f,(n) 1s approximated as f,.,(n) is 
expressed as | 

Ax+) (Nn) =f, (n)-2% w(n — 2m) fxs; (m). (4.67) 
The series expression on the right-hand side of (4.67) is recognizable as the same form of 
equation as appeared in (4.57). The form represents expansion of a sequence followed by FIR 
filtering. The factor of two is introduced to cancel the factor of 1/2 introduced by the 
decimation operation as indicated by (3.6). As with Mallat’'s algorithm, iteratively repeating, 


results in the cascaded tree structure of Figure 4.7a. Reconstruction of a decomposed signal 1s 
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Figure 4.7a--Block diagram illustrating three-stage decomposition by Laplacian pyramid. 
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Figure 4.7b--Block diagram of structure for reconstructing sequence decomposed by three-level decomposition 
depicted by block diagram of Figure 4.7. 
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Figure 4.8--Lattice illustrating concept behind generation of kernel w(n) used for Laplacian pyramid 
decomposition and reconstruction. After (30). 
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accomplished by simply reversing (4.67) and reintroducing to f,.,(n) the detail A, (n) removed 


during decomposition: 

f,.(n) = Axsi(n)+2-Zw(n—2-m) fy. mm). (4.68) 
Figure 4.7b illustrates the process by which a sequence decomposed by Figure 4.7a is 
reconstructed. 

Figure 4.8 provides a graphical illustration of the averaging process employed in the 
Laplacian pyramid. The Laplacian pyramid decomposition involves the approximation of some 
sequence f,(n) by an averaged version f,(n). The technique developed by Burt and Adelson 
[30], in calculating the values of the nodes in Pca endeavors to consider each node in f,(n) with 
an equal weight. In the case of a five-point weighting sequence w(n), each node in f,(n) is 
calculated from an average of its five nearest neighbors in f,(n). For example, the value of node 
f,(n) in Figure 4.18 is evaluated as 

f,(0) = w(-2)-fy(-2) + w(-1)-fy(-1) + w(0)-fy(0) + w(1)-fy( 1) + w(2)-f(2). 
Additionally, obtaining f,(-1) and f,(1), requires evaluation of 


f,(-1) = w(-2)-fo(-4) + w(-1)-f(-3) + w(0)-fo(-2) + w(1)-fo(-1) + w(2)-fo(0) 
and 


f,(1) = w(-2)-f5(0) + w(-1)-fo(1) + w(0)-fo(2) + w(1)-f9(3) + w(2)-f9(4). 
Now f,(0), an even-numbered node, contributes to the computation of three nodes in f,(n): to 


f,(0), to which it is directly adjacent across scale, and to f,(-1) and f,(1). However, f)(1), an 
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odd-numbered node, only contributes to the contribution of two-nodes in f,(n): to f,(0) and to 
er). 


As stated previously, the objective is to define {w(n)}_ such that the total weighted 
contribution by each node in f,(n) to the computation of nodes in f,(n) is equal. Therefore, the 
requirement that the total weighting factor for all contributions by f,(0) to equals the total 
weighting factor for all contributions by f,(1) implies that 

w(-1) + w(1) = w(0) + 2-w(2). (4.69a) 
The left-hand side of (4.69a) represents the sum of the weighting factors for all contributions by 
f,(1) to all nodes to which it contributes in f,(n). Similarly, the right-hand side of (4.69a) 
indicates the sum of the weighting factors for all contributions by f,(0) to all nodes to which it 
contributes in f,(n). Inductively, for a five-point weighting sequence the sum of the weights for 
all contributions by any even-numbered node in f,(n) will equal the right-hand side of (4.69a) 
while the sum of the weights for all contributions by any odd-numbered node in f,(n) will equal 
the left-hand side of (4.69a). 

An additional constraint was imposed on the selection of the weighting sequence w(n): 

= wk) = ly (4.69b) 
The constraint imposed by (4.69b) ensures that, for any node in f,(n), the sum of the weights by 
all contributing nodes in f,(n) will be unity. Consequently, the average energy in f,(n) 1s 


preserved by the approximation f,(n). 
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To satisfy (4.69a) and (4.69b), Burt and Adelson selected 


w= le 
i 4 as 


ie (4.70) 


tu |e 


Ls 
4 4 
Each odd-numbered node contributes twice, each time with a weighting of 1/4. The 
even-numbered nodes contribute with a weighting of "a" when directly adjacent across scale 
while contributing twice with a weighting of 1/4 - a/2 when separated by scale and translation. 


Consequently, each node contributes with a total weight of 1/2. 


Considering the top branch of the structure in Figure 4.6 from a two-sided Z-transform 
perspective provides insight into the frequency-domain character of the Laplacian pyramid. The 
transform of the weighting sequence can be expressed as 


Wiz) = (4-3) 42%) +h etz") +a. (4.71) 


Applying (3.6), decimating (4.71) produces 


Wi(z) = >[W@!?)+W(-z!?)] 


(L-2)@+24)4 a 


which after expansion becomes 
[Wan = (4-4) +27)+ a. (4.72) 
The effective transfer function T(z) of the top branch of the structure in Figure 4.6 can be 
expressed as 
T(z) = 2:[W,,],.(z)-W(z). (4.73) 
Through trial and error, Burt and Adelson observed that for the choice of a = 3/10, repeated 


applications of the transfer function (4.73) produced an impulse response with approximately the 
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same shape with each iteration and which ultimately, approached a Gaussian shape. For this 


choice of "a", the system transfer function becomes 


+ 28+4+27+42° + Ce + 2) ez +t 


T(z) = cee ee ee 


z3 


Furthermore, from Figure 4.6, by inspection, the expression for the detail sequence A, ,,(n) is, in 


terms of f,(n) and t(n), 


Acei(n) = % t(n —k) fe(k) — f(n). 
Equivalently, if the transfer function for the entire structure of Figure 4.6 is defined such that its 
impulse response is d(n), then d(n) = t(n) - 6(n) where 8(n) is Kronecker's delta function. 


Consequently, the transfer function D(z) becomes 


as a! ey eee ee ee eee. tk | 
-_— += -_— —_— 4 +—= + — 7+— 
Re zZ 5 z 5 zZ 0 Z 52 5 a 2 


Diz3 = ——— XA 


z4 


From this analysis follows a conceptually clearer but computationally less efficient structural 


equivalent for Figure 4.6. This equivalent structure is indicated in Figure 4.9. 





Figure 4.9--Equivalent structure to that shown in Figure 4.6. 


As a method of time-frequency decomposition, the Laplacian pyramid provides 
performance for processes containing high-frequency components which inferior to that which 


shall later be observed for Mallat's algorithm. Figure 4.10 illustrates the partitioning of the 
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frequency spectrum which occurs when a signal is decomposed by the Laplacian pyramid. The 


spectral bins which are formed do not follow the pattern of an even division of the frequency 
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Figure 4.10--Partioning of frequency spectrum perforned by three-stage Laplacian pyramid of Figure 4.7. 


spectrum at each stage as would be anticipated by operations involving factor-of-two decimation 
and expansion. In fact, approximately three fifths of the frequency spectrum below the Nyquist 
frequency is contained within the highien frequency'bin. Consequently, for classification 
applications, a Laplacian pyramid-based analyzer would provide poor localization of frequencies 
above approximately 0.2 f,. However, for an over-sampled sequence comprised primarily of 
lower frequency components, a Laplacian pyramid detector would likely provide slightly better 
spectral resolution than would be afforded by a multiresolution detector which divides the 
spectrum evenly at each stage. For a three-stage, even-division multiresolution scheme, the 
region in the frequency domain from 0 to 0.125-f, would be divided into two spectral bins. For 
the scheme of Figure 4.10, the region in the frequency domain from 0 to approximately 0.15-f, is 


divided into three spectral bins. 
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Figure 4.1 1--Block diagram for full expansion by channel approach for reconstruction of a sequence decomposed 
by the Laplacian pyramid. Structure is equivalent to that of Figure 4.7b. 


In general, multiresolution decomposition schemes which involve factor-of-two 
decimation generate lattices which resemble Figure 4.8. At each stage, the approximation 
density is reduced by a factor of two. For information transmission systems, this represent the 
primary advantage of such systems. However, for display of time-scale decompositions 
involving mesh or contour plots, most computer graphical routines require the insertion of zeros 
in order to create a lattice of points of constant density. For mesh plots, such as those provided 
by Matlab, a lattice of the form of Figure 4.8 will appear as rows of fin-like structures of 
varying density. Because of the rapid transition to zero at each lattice nodes, contour plots 
essentially provide binary indications of a sequence’s time-scale content: Indication of the 
presence of energy at a node is given, without indication of the relative magnitude of the 


sequence at that node. 


Figure 4.11 provides a block diagram of a full expansion by channel scheme for 


preparing the lattice similar to that of Figure 4.8 for plotting. The structure of Figure 4.11] 


87 


represents the result of separating each branch of Figure 4.7b at its summation points. The 
consequence of implementing a structure like that in Figure 4.11 is that sparse rows in a lattice 


like that in Figure 4.8 will be interpolated resulting in a uniform density of lattice nodes. 


Sequence value, s(n) 
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Figure 4.12--Time plot of 256-point test sequence generated by (4.74a) for demonstration of Laplacian Pyramid. 
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Figure 4,1 3--Plot of power spectral density of 256-point test sequence generated by (4.74a). 


Furthermore, the reconstructed sequence from the structure in Figure 4.7b is 


accomplished by 
fo(n) = x AM) (1), 
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In other words, the reconstructed sequence fo(n) consists of the sum of the content of all the 
channels f*™(n) obtained from Figure 4.11. 


To demonstrate the operation of the Laplacian pyramid, two test sequences were 


employed. The first sequence was constructed from low-frequency components: 


SiF(nN) = (1 | enya. e275 | | 2 -cos(2-m >n)+ sin(2-m-S-n)++ cos(2-1- sin)]. 


(4.742) 
A time plot of the sequence is presented in Figure 4.12 and a power spectral density plot appears 
in Figure 4.13. The test sequence was chosen such that the spectral content would lie 
predominantly in the lower portions of the frequency spectrum. Test sequence (4.74a) contained 


one spectral component, which completes only one complete oscillation during the duration of 


the sequence, at 7/128 (= 0.004-f,). Additional components are present at 9:7/256 (= 0.018-f,) 


and at 31:2/256 (= 0.06'f,). 


The second sequence, identical to (3.47) employed in Chapter III, Section C 1s repeated 


here for convenience: 


—{ |] —~ e764 |. o-{n-94)21 , 920/75 | eae ll ae. li ee lay ] 
sur(n) = (1 e e e [ cos(2 Mase) Meh cos(2- 7 ee Hig cos(2- 7 ea Tl) 


(4.74b) 
A time plot of (4.74b) is presented in Figure 3.23 and its power spectral density is plotted in 


Figure 3.24. Test sequence (4.74b) contains spectral components at 2:%:27/256 (= 0.105-f,), 


2-m-55-1/256 (= 0.215-f.) and at 2:2-113-7/256 (= 0.44-f,). 
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Figure 4.14a presents a surface plot depicting the eight-scale, full expansion by channel 
of the Laplacian pyramid decomposition of the test sequence generated by (4.74a). Despite the 
fact that the spectral components of the lowpass test sequence (4.74a) are quite close together, 
the Nevertheless, the Laplacian pyramid 1s able to distinguish between the three separate 
components. One peak 1s evident at around sample 32 at a scale of five to six, a series of 
approximately four peaks are evident at a scale of three to four from samples 10 to 60 and a 
series of peaks corresponding to the high-frequency components are apparent at scales one to 


two from samples 20 to 120. 
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Figure 4.14a--Surface plot of full expansion by channel of Laplacian pyramid decomposition of 256-point, 
lowpass test sequence generated by (4.74a). 


Figure 4.14b illustrates the Laplacian pyramid decomposition of the highpass 256-point 


test sequence (4.74b). As discussed within the context of the frequency partition diagram, 
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Figure 4.10, the Laplacian pyramid produces significantly poorer resolution for sequences 
comprised of high-frequency components. In fact, for the highpass test sequence (4.74b), the 
two spectral components with the highest frequencies fall within a single spectral bin. 
Consequently, as Figure 4.14b indicates, it is very difficult using the Laplacian pyramid to 


resolve spectral components residing in the upper two thirds of the frequency spectrum below 


the Nyquist frequency. 
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Figure 4.14b--Surface plot of full expansion by channel of Laplacian pyramid decomposition of 256-point, 
highpass test sequence generated by (4.74b). 


The reconstruction error for the Laplacian pyramid operation on the test sequence was 
consistently outstanding. The mean square reconstruction error as defined by (3.32) is plotted in 


Figure 4.15. As indicated by Figure 4.15, the reconstruction error never rose above -317 dB. 
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This accuracy roughly corresponds to the numerical precision of Matlab which was used for this 


demonstration. 


Mean Square Reconstruction 
Error (dB) 





2 3 4 5 6 7 g 
Number af Oecomoasition Staaes 


Figure 4.15--Trend of mean square reconstruction error for Laplacian pyramid operating on 256-point lowpass 
test sequence of (4.74a). 


The a trous algorithm represents the first multiresolution analysis technique explicitly 
based on an affine-type representation vector. This technique consists of evaluating a 
discretized approximation of the continuous wavelet transform. Specifically, the starting point 
for the a trous method Is the discrete wavelet series [31]: 


w(m, n) =27™ zwe(2™ k —n) s(k). (Ais 


Essentially, the discrete wavelet series consists of the discrete correlation between some 
sequence s(n) and translations of a sampled wavelet function 2™ w(2'k - n) at scale m. 


Evaluating (4.75) for w(1, n) produces 


w(1, n) = > Ly'(F-n) s(k). (4.76) 


In (4.76), it is obvious that w(t —n)= Td Consequently, (4.76) becomes 


w(1, n) = Dy" (S*) s(k). 
2 ok 


; 


v 
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Now, if y(k) is known for integer-valued arguments, its values at half-integer arguments 
can be approximated through interpolation. If f(k) represents the impulse response for an 


interpolation filter, then y(k) at half-integers is approximately 


w(4) = /2 ZF (2k-n) y(k), (4.77) 


where the conjugation and time-reversal of the filter and the factor of /2 have been introduced 
for later convenience. The form of the summation term on the right-hand side of (4.77) is 
recognizable as expansion of a sequence w(k) followed by convolution with a filter whose 
impulse response is f(-k). 

To obtain the best approximation, special conditions are imposed on the filter represented 


by f(k). Shensa [31] applies the term a trous filter and employs the notation 


f{(2-k)= = do, 


where 5, , is Kronecker's delta function. Perhaps a slightly more illuminating manner in which 


to express this concept is to use the notation 


fy.(k) = F Bok: (4.78) 
In other words, the decimated impulse response f(k) 1s non-zero only fork =0. Equivalently, 
the only non-zero, even-numbered coeffecient for the filter f(k) corresponds to k =0. The 
condition (4.78) is necessary in order to ensure that known values for the interpolated sequence 


w(k) are unaffected by the interpolation process. When evaluating f (2k - n) w(k), ifn is an 


even integer, an argument for which y(n/2) is known, 2k - n will also be an even integer, and 


a3 


the series (4.77) will contain only one non-zero term, for n= 2k. Consequently, (4.77) will be 


evaluated as 


Hence, the factor of ae in (4.78). On the other hand, if n is an odd integer, an argument 


value for which the interpolation must be calculated, the series (4.77) will be evaluated as a 
weighted sum of the function values for y(n/2) for surrounding even values of n. The a trous 
condition (4.78) assumes a position of importance tn the theory of FIR filters. The condition 
(4.78) is equivalent to the half-band condition. A Aalf-band filter is a symmetric, FIR digital 
filter whose impulse response satisfies (4.78). “Mitzner [32] showed that filters satisfying (4.78), 
of necessity, also satisfy 

F(e!®) +F(e!!") =2f(0) where F(e!) = ~ fig”. (4.79) 
Consequently, the interpolation filter whose impulse response ts f(k) is also a half-band filter. 
Furthermore, Shensa showed that the filters for Daubechies’ orthonormal scaling functions 
belong to a special class of filters satisfying (4.78) known as Lagrangian interpolation filters. 
Specifically, if c(n) 1s the impulse response for the filter producing one of Daubechies' 


orthonormal scaling functions, the corresponding Lagrangian interpolation filter 1s obtained from 





f(n) = = E c(k) e(k eit 


Continuing the development of the a trous algorithm, substitution of (4.77) into (4.76) 


produces 


w(1, n) = Xd f(2j+2n—k) w°()) s(k). (4.80) 


j} k 
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Now, if a filter g(n) 1s defined such that 
g(n)=ay-(-n), (4.81) 


and a change of indices m =j + nis applied to (4.80), it follows that 


w(l, n) =2g(n—m) 2 f(2m—k) s(k). (4.82) 
Equation (4.82) simply describes decimation of s(k) followed by filtering with a FIR filter 
whose impulse response is f(k). The result is then further filtered by an additional FIR filter 
whose impulse response is g(n). Furthermore, the second FIR filter is the conjugated 


time-reversal of the sampled wavelet function. Therefore, filtering with g(n) is equivalent to 


evaluating the correlation with y(n). 
- | 


fi in) 


Figure 4.16--Block diagram of one-stage decomposition by a trous algorithm. 






Now, the operation of decimation and then filtering s(n) resembles the portion of Mallat's 
algonthm, (4.53), where successively coarser approximations of s(n) were obtained by 
decimating and filtering each preceding approximation. For the signal detail is extracted merely 
by filtering the upper half of the spectrum below the Nyquist frequency. If the notation s'"(n) is 
adopted for the i*-level approximation of s(n) and w(n) for w(i, n), by induction, from (4.82) it 
follows that 


Se ee, (4.83) 
w")(n) =g * s(n) 
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The block diagram for (4.83) 1s shown in Figure 4.16. 


Earlier wavelet transform techniques resemble Gabor transforms in which sequences 
were decomposed by projection on representation vectors consisting of modulated Gaussian 
windows. Consequently, the wavelet function used by Shensa in his study of wavelets was a 
modulated Gaussian, also known as a Mor/let window: 

y(t) =e Vee Fi? (4.84a) 

Fourier transformation of (4.84a) produces 
Wo) = 7 2m eee (4.84b) 
In order to specify a wavelet function, it is, therefore, necessary to specify its center frequency v 


and its window rolloff factor B. 


Furthermore, in order to increase resolution in the Fourier domain it is possible to 
employ a wavelet function consisting of superimposed Gaussian envelopes modulated at 
different frequencies. In this case, the wavelet y(t) and its Fourier transform become 


M-1 (4) erin 
w(t) = ry e Z ako e sho 
SS see (4.85) 
—-Vj)° 
M 


, et k= 
(w)=4 /2m XD 2Me 2B 
B 
k=) 





Multiple frequency-domain translations of the Gaussian window are referred to as voices. 
Through the use of voices, the upper half of the frequency spectrum below the Nyquist 
frequency 1s partitioned into multiple spectral bins. Shensa listed a series of guidelines for 


designing a wavelet function with multiple voices: 
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The center frequency must lie in the upper half of the frequency spectrum below the 
Nyquist frequency. Consequently, 
n/2 <v. (4.86) 
In order for y(t) to be an analytic function, 
—_. (4.87) 
Obviously, (4.84b) has an infinite region of support. However, if, for (4.87), the equality 


is selected, 
Ww(0) = /2n e?™ =6.71x 10>. 


Consequently, the spectral content in regions of negative frequency will be negligible. 


To exclude aliasing, 


vsrn—/2 8B. (4.88) 


The quantity /2 Bis equal to one half of the Gaussian window where, at the edge of the 
passband, the frequency response of the filter is less than its maximum by a factor of e”. 
Condition (4.88) ensures that the high-frequency edge of the passband is located below 


the Nyquist frequency. 


The fourth condition for designing wavelets employs the concept of relative bandwidth. 
The relative bandwidth B, of the window is defined to be the ratio of the window 


bandwidth to its center frequency, or 





Since, from (4.86) and (4.88), m/2 <v <7, it follows that 


2 


ene AEE Bi (4.89) 
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tan 


This condition is essentially a combination of previous conditions. 


Furthermore, if multiple voices are desired, it follows that the number of voices M 


required to completely cover the upper half of the spectrum is 


tw 


M=— 


¥ 





— 
pe 2 Pe 


tv | 


tw 


As a result, an approximate choice for B is 





— | 
ps (4.90) 
To construct the separate voices, the mother wavelet y(t) is scaled by various factors 2*™ 


where M indicates the number voices. Specifically, the k” voice v,(t) fork =0, ..., M-1 


is given by 
V,(t) =v(str] 4.9 1a) 
from which it follows by the scaling property of the Fourier transform, 
¥4(@) = F2™ 20 eM o-uit2 BF (4.91b) 
This scaling produces the combined affects of increasing the Gaussian window decay 
factor by a multiple of 2*”™ and of translating the location of the window's maximum 


value to @ v. Selection of this scheme also maintains the affine nature of the 


transform. 


Then, to select the filter length, it is appropriate to consider the decay of the Gaussian 


window in the time domain. The filter will consist of a discretely sampled version of the 


continuous wavelet function. The DFT of the wavelet is, by definition, 


Vitel = 2. v.(k) e7JO™, 
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Due to the rapid decay of the window function y(t), the series can be truncated for values 
of the summation index k such that v,, ,(k), the voice with the slowest time-domain rate 
of decay, will be negligible. [fa negligibility threshold T, 1s specified, then a filter 
half-length K must be selected such that 


en ulead 


Inverting the envelope portion of (4.91a) produces for the filter half-length K 


K > ne Cie /—) in(T,) . 92) 


Consequently, the analyzing wavelet filter g(n) from (4.81) is defined to be 


*(- Vv a ce CK 
ee a (4.93a) 
0 otherwise 


and the filter g,(n) corresponding to the k" voice v,(t) becomes 


g.(n) = w*(—n/2EM) VY n=-K,..., K 
3 0 otherwise 
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Figure 4.1 7--Spectral partitioning perfonned by three-stage decomposition by a trous method using four-voice 
Morlet window as analyzing wavelet. The interpolation filter, plotted as the lowpass filter represents a 
Lagrangian filter calculated from Daubechies' orthononnal scaling function on (0, 3]. 
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To illustrate the spectral behavior of the a trous method, Figure 4.17 plots the frequency 
decomposition performed by a three-stage structure employing a four-voice Morlet window. 
The spectral bin density increases in a logarithmic manner as digital frequency approaches zero 
and the width of each bin decreases similarly. The relatively large bandwidth of the lowest 
trequency bin occurs because of the relatively small number of stages. Additional stages would 
result in the continued trend of logarithmically increasing spectral bin density and bin width. 
The logarithmic trend must, however, be ultimately broken by the final spectral bin which will 


be defined by a half-band, lowpass filter. 


Figures 4.18a and 4.18b demonstrate the performance of the a trous algorithm for 
identification of the spectral components comprising a sequence. Figure 4.18a presents the 
decomposition of the lowpass test sequence (4.74a). Examining Figure 4.18a, the three 
harmonic components can be observed. Ata scale of approximately 2.5 to 3, the band 
corresponding to the highest frequency component is evident. The middle frequency component 
1S apparent in the range from three to four. The lowest frequency constituent ts indicated by the 
sparce group of fins at scales of five to seven. The shape of the envelope is evident from the 
relative heights of the peaks of the two components possessing higher frequencies. However the 
lower frequency component appears relatively constant across the entire range of samples. This 


occurs because the integration time of the final scale extends the duration of the sequence. 


The a trous decomposition of the highpass test sequence (4.74b) ts plotted in Figure 
4.18b. Unlike the Laplacian pyramid, the a trous method displays no difficulty in resolving the 


spectral components of (4.74b). Figure 4.18 provides an opportunity to understand the 


1QQ 


relationship between the location of spectral components in the time-scale plane and their 


location in the frequency spectrum. The regions contained within each spectral bin are, again, 


indicated in Figure 4.17. The spectral component of the test sequence at 0.44-f falls into the 


right-most spectral bin of Figure 4.17. This component 1s indicated with in Figure 4.18b at a 


scale of zero from approximately samples zero to 160. The location of the spectral component at 


0.21-f, approximately coincides with the location of the boundary of the fifth and six spectral 


bins from the right-hand side of Figure 4.17. Since the analyzing wavelet contains four voices 
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Figure 4./8a--Surface plot of a trous wavelet transform decomposition of 256-point, lowpass test sequence 
wenerated by (4.74a). Plot represents eight-stage decomposition using four-voice Morlet analyzing wavelet and 
Lagrangian interpolation filter calculated from filter for Daubechies' orthonomnal scaling function supported on 
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Figure 4.18b--Surface plot of a trous wavelet transform decomposition of 256-point, highpass test sequence 
venerated by (4.74b). Plot represents eight-stage decomposition using four-voice Morlet analyzing wavelet and 
Lagrangian interpolation filter calculated trom filter for Daubechies' orthonormal scaling function supported in 
[0.0312 


(therefore there are four spectral bins per scale), the fifth and sixth spectral bins coincide with 
scales of 1.0 and 1.25, respectively. These scales coincide with an indication on 4.186 from 
samples zero to 96. Finally, the spectral component at 0.105-f, is located in the 10" spectral bin 
of Figure 4.17. This bin coincides with a scale of 2.25. Therefore, this final component, 
possessing the greatest power, provides an indication in Figure 4.18b throughout the entire range 


of samples. 


In working with the a trous algorithm, one significant disadvantage arises. Specifically, 


reconstruction of a sequence decomposed by the a trous method is extremely difficult. The 
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wavelet basis functions used in the a trous decomposition comprise frames in the sense described 
in Chapter II. Consequently, in order to reconstruct a sequence decomposed by the a trous 
method, it is necessary to construct biorthogonal set of bases. Daubechies outlined the details of 


construction the dual of a specified basis in [6]. 


Apart from the difficult of inversion, the a trous algorithm, as implemented by Shensa, 


suffers from the disadvantage of long filter lengths. In the example under consideration, in order 


to satisfy a negligibility threshold (4.92) of 5 x 10°, a filter length of 81 was necessary. Each of 
the four voices was implemented with a filter of equal length. Furthermore, the wavelet filter 
coefficients are complex. Complex filter coefficients increase the computational burden. 
Furthermore, from Figure 4.17 it is evident that the spectral decomposition resulting from the-a 
trous method is anything but power complementary. Partitioning in a manner similar to that of 
Figure 4.17 results in exaggeration of spectral content near the peaks of the voices and subdued 
indications between the peaks. In favor of the a trous algorithm, however, 1s the attainability of 
arbitrarily high resolution with smooth, regular filters. (It will be observed in the next section 


that filters sometimes suffer from lack of regularity.) 


E. MULTIRESOLUTION ALGORITHMS FROM CASCADED FILTER BANKS 


Section C of this chapter outlined the equivalence between Mallat's multiresolution 
analysis and a structure constructed from cascaded, two-channel, perfect reconstruction QMF 
banks. This section explores, in a more general manner, implementation of multiresolution 
analysis structures constructed from cascades of filter banks of arbitrary numbers of channels. 


For each method, the issues of division of the frequency spectrum and invertibility will be 
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considered. Spectral properties of the structures are of interest to obtain an understanding of the 
frequency resolution attainable with techniques under consideration. Invertibility is important 
because, in general, one good test of a decomposition technique is the accuracy of its 


reconstruction. 


a,(n) b(n) 
yee b(n) 





Figure 4.19--Block diagram of three-stage, two-channel multiresolution structure constructed from cascaded QMF 
analysis banks. Algorithm represented 1s equivalent to Mallat's algonthm (4.53). 


A block diagram of the most basic multiresolution algomthm is depicted in Figure 4.19. 
The process of Figure 4.19 1s exactly equivalent to that of (4.53), Mallat’s multiresolution 
decomposition. At each stage, the approximation sequence a,(n) is filtered with high-pass and 
low-pass half-band filters. The output of each filter 1s then dilated. The low-frequency channel 
decimator output--or the approximation channel, conforming to the noe of 
multiresolution theory is then applied to another stage. The high-frequency, or detail, channel is 
employed for transmission, in the case of a communication system, or is applied to a detector in 


one possible case of a digital signal processing application. 
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If, hypothetically, the filter for each channel in the structure depicted in Figure 4.9 is an 


ideal filter, the DFT filter output sequence for the approximation sequence bandlimited to 
L*({-n/2,n/2]), applying concepts from Chapter III. The DFT of the detail sequence is likewise 
bandlimited to L*((/2, 3-1/2]). Decimation of each sequence dilates its spectrum so that the 


approximation and detail sequence spectra are members of Caer m]) and L?([r, 2:1}), 


respectively. Consequently, the time-domain decimation and resultant frequency-domain 


dilation result, at each stage, in the examination of a smaller segment of the frequency spectrum. 
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Figure 4-20--Partitioning of the frequency spectrum resulting from five-stage decomposition structure of Figure 
4.19. Structure was implemented using, for the approximation and detail channels, the filters corresponding to 
Daubechies' orthonormal scaling function and wavelet, respectively, supported on [0, 13]. 


Figure 4.20 plots the spectral partitions resulting from a five-stage decomposition 
structured as in Figure 4.19. The spectral divisions indicated represent those obtained using the 
filters for Daubechies' orthonormal scaling function and wavelet on [0, 13]. The logarithmic 
contraction of succeeding spectral bins constitute the multiresolution analysis. The spectral 


partitioning of Figure 4.20 compares favorably to that obtained by the Laplacian pyramid in 
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Figure 4.10. However, the localization provided by the partitioning of the frequency spectrum 
plotted in Figure 4.17 obtained by the a trous is obviously vastly superior to that of 4.20. For 
some applications, the resolution of Figure 4.20 may be inadequate while the computational 


burden necessary to obtain the results of Figure 4.17 is excessive. 


Brooks [33], in experimentation with signal processing techniques for ultra-wideband 
radar systems, observed that a dichotomous spectral decomposition, as is performed by the 
structure of Figure 4.19, provided poorer results than the a trous algorithm. One possible reason 
for the disparity between the performance of the two techniques considered was the inadequacy 
of the spectral resolution provided by Mallat's algorithm. Consequently, consideration of a way 
in which to improve the spectral resolution of multiresolution techniques without incurring the 


disadvantages of the a trous algorithm 1s justified. 
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Figure 4.2/--Multiresolution structure obtained by cascading detail channel of Figure 4.20 with a three-channel 
QMF-type filter bank. 












Figure 4.21 represents one possible approach to improving the results of Figure 4.20. 
Applying, at each stage, the detail channel to the input of the analysis bank of another filter bank 


produces a signal decomposition method capable of obtaining arbitrary precision. The study of 
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filter bank design is an advanced field and techniques for the design of perfect reconstruction 
filter banks of an arbitrary number of channels are available employing lattice structures [11] or 


through cosine modulation of a properly designed prototype filter [34]. 


Figure 4.21 presents one example of a single stage of a structure using subchannel 
decomposition of the detail channel at each stage. In Figure 4.21, at each stage, the spectral bin 
corresponding to the detail sequence is further subdivided into three subchannels. Figure 4.22 
illustrates the spectral partitioning which results from a four-stage structure constructed from 
stages indicated by Figure 4.21. The improvement in spectral resolution is obvious. The density 
of groups of subchannels increases as digital frequency approaches zero. However, the overall 
structure does not provide the logarithmic spacing of spectral bins observed for the a trous 


decomposition as indicated in Figure 4.17. 


All 


hy soaks NNN OK 


Wiad lac 


A\s AS WAY \ i Px) ites 


@ 2.05 @.1 @515 @.2 00025 Od 2.959 G4 O45 @9 
Digital frequency (multiple of tied 


——=— 


ol Ov Co 


Normalized Magnitude 


NO 


Figure 4.22--Partitioning of frequency spectrum resulting from four-stage multiresolution structure constructed 
from stages appearing in Figure 4.22. Structure was implemented using, for the approximation and detail 
channels, the filters corresponding to Daubechies’ orthonormal scaling function and wavelet, respectively, 

supported on [0, 13]. Subchannel filter banks were implemented with three-channel filter bank coefficients are 

presented in [11]. 
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Subchannel decomposition can be performed using either perfect reconstruction QMF 
banks or pseudo-QMF banks--filter banks which do not strictly satisfy the perfect reconstruction 
criteria. Similarly, to the extent in which the pseudo-QMF bank deviates from the power 


complementary property, accurate representation in scale will be forfeited. 





Figure 4.23--Structure constructed from cascade of three-channel filter banks. 


An additional method to improve spectral resolution of multiresolution techniques was 
suggested by Gopinath and Burrus [35]. Employing the term mu/tiplicitvy M wavelet transform, 
Gopinath and Barrus suggested the concept illustrated in Figure 4.23 as a generalization of the 
dichotomous structure proposed by Mallat. In the frequency domain, spectral resolution 
Increases as a base three logarithm instead of the base two logarithmic contraction observed in 
frequency spectra resulting from Mallat-type structures. Furthermore, structures similar to those 


of Figure 4.3 can be cascaded from filter banks of any number of channels. 
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Figure 4. 24--Spectral partitioning resulting from four-stage decomposition constructed of structures as in Figure 
4.23. The filters indicated are 14-point pseudo-QMF filters designed in [11]. 


Figure 4.24 illustrates the partitioning of the frequency spectrum which results from a 


b) 





Figure 4.25--{llustration of equivalence of two-stage decomposition-reconstruction structure with one-stage 
structure with delays in each channel. 


In reconstructing sequences decomposed by cascaded filter banks, the issue of delay must 
be addressed. As was discussed within the context of multirate system theory in Chapter II, a 
perfect reconstruction filter bank possesses an equivalent system transfer function which 1s equal 
to adelay. In decomposition-reconstruction structure "a" of Figure 4.25, the detail channel 


contains a delay L, representing, for example, a reconstruction filter bank used, as in Figure 
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4.21, to decompose the detail sequence into subchannels. The lower branch, the approximation 
channel, contains a two-channel bank corresponding to the next stage of decomposition. In the 
equivalent structure in structure "b” of Figure 4.25, the follow-on stage filter bank has been 
converted to the delay L, with which it is equivalent if it 1s a perfect reconstruction filter bank. 
Obviously, the delay in each channel must be equal. Otherwise, the approximation and detail 
sequences will be shifted with respect to each other and accurate reconstruction will not occur. 
Consequently, although a subchannel decomposition filter bank 1s not indicated in the detail 
channel of the final decomposition stage in Figure 4.25, if one were installed, an artificial delay 


in the approximation channel would be necessary. 


zbs 
a) 


) zbs 


b) 





C) 


Figure 4.26--Ilustration of equivalent structures resulting from transmission of delay occurring in approximation 
and detail channels of a filter bank mulstiresolution structure. 


Figure 4.26 illustrates the total equivalent delay of a QMF decomposition-reconstructure 


containing a delay in each channel. In each channel a delay L, (equal to the order of the filters 
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used for decomposition into subchannels) is indicated which could include combinations of 
delays of subchannel decomposition-reconstruction structures, delays from subsequent 


decomposition-reconstruction stages and artificial delays introduced to resynchronize channel 


feauences. Since, by (3.9), [z as =z°*'* the total effective delay after the expander in each 
channel will be twice the delay of the channel prior to the expander. This equivalence is 
indicated by the transition from structure "a" of Figure 4.26 to structure "b." Furthermore, since 
the application of FIR filter and application of a delay represent linear operations, the two 
operations are also commutative. Consequently, as indicated in structure "c" of Figure 4.26b, 
the total, effective delay of the channel after the decimator can be moved beyond the summation 
point. Since a sequence reconstructed by a perfect reconstruction QMF bank is simply a delayed 
version of the original sequence, the total equivalent delay of structure "a" in Figure 4.26 is 


simply the sum of L, (the order of the filters used for decomposition into primary channels), the 


delay of the decomposition-reconstruction QMF bank structure itself, and 2-L,, the contribution 


from all of the delays introduced into the decomposed channels. 


To obtain a characterization of the total delay introduced into each channel of a 
multiresolution analysis constructed from cascaded filter banks, consider again structure "a" in 
Figure 4.25 as it would be implemented to reconstruct a sequence decomposed by M stages of 
the structure depicted in Figure 4.21. First, if a delay of L., resulting from a subchannel 
decomposition-reconstruction structure, is introduced into the detail channel at stage M, of the 


final stage of decomposition, then an artificial delay of L, must be introduced into the 


approximation channel of the same stage to resynchronize the detail and approximation channel 
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sequences. Consequently, the total delay introduced into the detail and approximation stages 
consists of the sum of two contnbutions. The first source includes the delay L, resulting from 


decomposition-reconstruction structure corresponding to the decomposition of the detail 


sequence into subchannels and its corresponding equal delay in the approximation channel. 
Secondly, a delay L, corresponding to the primary channel decomposition-reconstructure at 


stage M must be considered. 


Continuing to the next stage, recombining the approximation and detail sequences from 
stage M-1 to form the stage M-2 approximation sequence adds a delay of L, to the doubled 
equivalent delay present into the approximation sequence of stage M-1. As a result, the total 
delay present in the approximation sequence at stage M - 2 is equalto4-L_+2:L.+L.. Itis 
therefore necessary to introduce a delay of 4-L, +L, + L,. to the detail sequence of stage M - 2 


to resynchronize that stages approximation and detail sequences. 
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Figure 4.27--Block diagram indicating structure necessary to obtain full expansion by channel from a sequence 
decomposed by M stages of the structure depicted in Figure 4.23. The values for L, and L, are defined in (4.94a) 
and (4.94b), respectively. 
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To obtain the structure required to obtain a full expansion by channel of the decomposed 
sequence by M stages of the structure presented in Figure 4.21, inductively the line of reasoning 


begun in Figures 4.25 and 4.26 must be pursued. The consequences of this induction appear in 
Figure 4.27. First, some arbitrary stage m, the subchannel sequence dee (a) must first be 


expanded and then filtered by the synthesis filter I,(z) corresponding to its subchannel. 
Artificial delays from two sources must then be introduced in order to re-synchronize the 


resultant stage m detail sequence component with its corresponding approximation sequence. 


The first delay L, compensates for the artificial delay of L, introduced into the approximation 


sequence at stage M and ts evaluated as 
M-m 
a 1 ES (4.94a) 


The second artificial delay contribution results from the primary approximation channel 
decomposition-reconstruction performed at all subsequent stages m+l through M and is 


represented by a geometric series: 


M-m-! 
be ccelee - ed, a22e (4.94b) 


p=0 
Next, the detail sequence component must be expanded and passed once through the 
primary detail channel synthesis filter. To complete the re-expansion, the sequence Is subjected 
to m-I iterations of expansion followed by filtering with the primary approximation channel 
synthesis filter. The final stage approximation sequence is fully expanded--as indicated by the 
second structure of Figure 4.27--by M iterations of expansion followed by filtering with the 


primary approximation channel synthesis filter. Application of the steps described by Figure 


i3 


4.27 provides for a decomposed sequence a representation analogous to that depicted in Figure 


4.11 for the Laplacian pyramid. The reconstructed sequence is simply the sum of all of the 


expanded channels or, if applicable, subchannels. 
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Figure 4.28a--Plots of reconstruction error versus number of decomposition stages for various decomposition 


stages applied to 256-point, lowpass test sequence produced by (4.74a). 


The accuracy of sequence reconstruction 1s dependent on the number of stages of 


decomposition performed and on the quality of the filter bank. For several stages of 


decomposition, small errors resulting from filter imperfections, such as roundoff error, 


accumulate and corrupt the reconstruction process. To obtain optimum results, filter banks 


which strictly satisfy the perfect reconstruction criteria should be used. In using pseudo-QMF 


filter banks which do not satisfy the perfect reconstruction, poorer results will be obtained. 


These concepts are illustrated by Figures 4.28a and 4.28b. Figure 4.28a indicates reconstruction 


eirors for multiresolution structures applied to the test sequence generated using (4.74a). 


The lowpass test sequence from (4.74a) was constructed from low-frequency harmonics. 


From the spectral partition diagrams presented for the various multiresolution structures 
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considered, the spectral resolution of the multiresolution structures was greatest in the lower 
regions of the frequency spectrum. In Figure 4.28b, the reconstruction error for the highpass test 
sequence (4.74b) is presented. Illustrating results for the test sequence generated by (4.74b), 
Figure 4.28b illustrates multiresolution technique performance for signals constructed from 
frequency components located in the upper half of the frequency spectrum below the Nyquist 


frequency. 
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Figure 4.28b--Plots of reconstruction error versus number of decomposition stages for various decomposition 
stages applied to 256-point test sequence produced by (4.74b). 


In the case of both test sequences, the lowest reconstruction errors resulted from 
multiresolution structures comprised strictly of perfect reconstruction filter banks. In 
two-channel, zero-subchannel and two-channel, two-subchannel structures, the filters used were 
all obtained from Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. 
Consequently, the multiresolution analyses obtained were equivalent to true, orthogonal 


decompositions. 
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For the two-channel, three-subchannel, the two-channel, four-subchannel and the 
three-channel, zero-subchannel cases, the structures all entailed the use of pseudo-QMF banks 
which did not strictly satisfy the perfect reconstruction criteria. The structures constructed using 
the three-channel filter banks employed the set of filters designed by Vaidyanathan in [11]. The 
four-channel filter bank used a set of 30-point, fourth-band filters designed by spectral 
factorization of filters obtained by the McClellan Parks algorithm. The design process for the 
four-channel filter bank will be outlined in Chapter V. Because of failure to satisfy the perfect 
reconstruction property, structures using these methods consistently produced greater 


reconstruction error. 


TABLE 4.1~PARTITION OF SPECTRUM RESULTING FROM EIGHT-STAGE, TWO-CHANNEL, - 
ZERO-SUBCHANNEL MULTIRESOLUTION DECOMPOSITION. 


Spectral bin Spectral bin 
lower bound Upper bound | 














Decomposition stage Scale 


Stage | Detail 





Stage 2 Detail 
Stage 3 Detail 
Stage 4 Detail 
Stage 5 Detail 
Stage 6 Detail 
Stage 7 Detail 
Stage 8 Detail 
Stage 8 Approximation 


Inductively extending the results of the spectral partition diagram (4.20) corresponding to 
this structure provides insight into the relationship between scale and frequency. The first 


spectral bin, which coincides with scale zero in Figure 4.29a, constrains the region of the 


frequency spectrum from 0.5-f, to 0.25-f,. A scale of unity coincides with the region of the 
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frequency spectrum from 0.25-f, to 0.125-f,, the region corresponding to the signal detail 
extracted during the first stage of decomposition.. Inductively, therefore, a scale of m, 
corresponding to the signal detail extracted at the m"" stage of decomposition, coincides with the 


spectral bin partitioning the region from ef to Per fs This subdivision is tabulated in 


Table 4.1. 
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Figure 4. 29a--Surface plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthononnal wavelet and scaling function supported on [0, 13]. 


Figures 4.29a and 4.29b demonstrate the full expansion by channel of decomposition of 
sequences (4.74a) and (4.74b), respectively. In Figure 4.29a, the depiction of the lowpass 


sequence (4.74a), the resolution of scale of the decomposition at low-frequencies is apparent. 


The component of sequence (4.74a) located at digital frequency 2--9/512 (= 0.018-f,) falls in 
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the range of 2 “ to 2° a which corresponds to a scale of four. In Figure 4.29a, a series of 
circular contours centered on the scale axis for m = 4 can be observed from approximately 
samples 32 to 216. Similarly, the spectral component located at digital frequency 2:2:31/512 


(= 0.06-f,) lies in the region 2 of to 2 -f. corresponding to a scale of m= 3. Asa non-integer 


power of two, the location 0.06°f, = De 


which 1s very near the boundary between adjacent 
scales two and three. Consequently, much of the component's energy appears spread between 


scales two and three. The component at digital frequency 2/256 (=0.002-f,) exists at the 


boundary between scales six and seven. The indication of this component, which completes 
only a single oscillation during the duration of the test sequence, appears spread between the two 


scales. 


Figure 4.29b presents the two-channel, zero-subchannel decomposition of highpass test 
sequence (4.74b). The relative coarseness of the multiresolution decomposition in the higher 
regions of the frequency spectrum becomes apparent in Figure 4.29b. Because of the rapid 
variations with respect to time, the contour plot in Figure 4.29b is difficult to read. However, 
information can be gained from examining the surface plot. Again, comparing the plot with the 
spectral partition diagram, the spectral component of (4.74b) which lies at a digital frequency of 
2:%:27/512 (= 0.105-f,) falls within the region of 2° * to 2°* corresponding to a scale of m= 2 as 


plotted in Figure 4.29b. The scale of m = 2 corresponds to the range of peaks with highest 


amplitude. The spectral component located at 2:2:55/256 (= 0.215-f.) is contained within the 


region 2°* to 2’* which corresponds to a scale of m = 1 as plotted in Figure 4.29b. The 
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indication of this component is recognizable in the second range of peaks between the highest at 


m = 2 and the range that runs along the edge of the surface plot. Finally, the spectral component 
at digital frequency 2:m-113/256 (= 0.44-f.) lies within the region of 2° ‘to 2°* and is indicated in 


the range of peaks at scale m = 0 which runs along the edge of the surface plot. For the highpass 


test sequence, each of the harmonic components can be resolved. However, in the case of the 


highest frequency components, the spectral separation is approximately 0.2-f,. Furthermore, 
each of the three spectral components fell within adjacent spectral bins. Consequently, the 


example demonstrates the maximum resolution for a two-channel, zero-subchannel 


multiresolution structure. 
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Figure 4.29b--Surface plot of decomposition of highpass test sequence (4.74b) using eighit-stage, two-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on (0, 13]. 
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Figure 4.30a--Surface plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
three-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 
performed with pseudo-QMF filters designed in [11]. 


Figures 4.30a and 4.30b provide two-channel, three-channel decompositions of the 


lowpass test sequence (4.74a) and the highpass test sequence (4.74b), respectively. From the 


palm.) saline*) 


spectral partition diagram for this structure, Figure 4.22, each spectral region to 
is further subdivided into three subregions. For the general case, a multiresolution structure 
constructed to divide, at each stage, the detail sequence into M subchannels, at the m" stage of 
decomposition, the k” subchannel corresponds to the region of spectrum contained in the interval 
| (2m + = ’ (220) = 2442) ata (2-00 ars ~ Cae ~ 2-42) | im 
Because of the increase in the density of spectral subdivision afforded by decomposition 


of the lowpass test sequence into subchannels, the time-scale localization presented in Figure 
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4.30a provides time-scale localization which is improved over results presented in Figure 4.29a. 
In particular, the component at digital frequency 2:2-31/512 clearly localized at a scale of three 
instead of being partially spread over adjacent scales as occurred in Figure 4.29a. Similarly, the 
component at digital frequency 2-7%:9/512, which roughly corresponds to a digital frequency of, 


af. can be recognized to exist primarily between scales of four and five. The component 
characterized by the lowest frequency 1s still recognizable in the vicinity of scales six and seven. 
However, an indication appears between scales of five and six which cannot be accounted for 
based on the components known to exist in the test sequence. This unaccounted for component 


may have appeared as a result of spreading from the component at scale six. 
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Figure 4.30b--Surface plot of decomposition of highpass test sequence (4.74b) using eight-stage, two-channel, 
three-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonomnal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 
performed with pseudo-QMF filters designed in [11]. 
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In the case of two-channel, three-subchannel decomposition, the spectral components of 
the highpass test sequence (4.74b) have also become more distinct. In the plot of Figure 4.30b, 
because the spectral content of the test sequence has been localized in specific subchannels and 
excluded from others, resulting in an indication of three separate components. Consequently, the 


indication provided by the contour plot constitutes a more useful display than the contour plot 


for Figure 4.29a. 
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Figure 4.31la--Surface plot of decomposition of lowpass test sequence (4.74a) using eight-stage, two-channel, 
four-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained from 
Daubechies' orthonormal wavelet and scaling function supported on [0, 13]. The subchannel decomposition was 

perfonned with pseudo-QMF filters obtained from spectral factorization of fourth-band filters designed by the 

McClellan-Parks method. 


Examining Figure 4.31a, a plot of a two-channel, four-subchannel decomposition of 
lowpass test sequence (4.74a), increasing structural complexity from three to four subchannels 


provides some improvement in localization for the test sequence presented. The component of 
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digital frequency 2:%-9/512 appearing between scales four and five have become distinct from 
the indication unaccounted for between scales five and six. Consequently, it is likely that the 
false indication resulted from spillover of energy from the component at 2:7:2/512. Between 
scales five and six, the resolution of the decomposition using four subchannels becomes 
0.25-2°F, (=f/512). Beneath those scales the resolution continues to increase logarithmically in 
base two. It appears, therefore, as though the four-subchannel decomposition may provide 
excessive resolution for the test sequence demonstrated. In general, the principal benefits of 


decomposing sequence into subchannels occurs in the higher regions of the frequency spectrum. 
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Figure 4.31b--Surface plot of decomposition of highpass test sequence (4.74b) using eight-stage, two-channel, 
four-subchannel multiresolution structure constructed from cascaded QMF banks. The filters were obtained trom 
Daubechies' orthonormal wavelet and scaling function supported on (0, 13]. The subchannel decomposition was 

performed with pseudo-QMF filters obtained from spectral factorization of fourth-band filters designed by the 

McClellan-Parks method. 
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Figure 4.31b further illustrates the principal afforded by subchannel decomposition. 
Specifically, in Figure 4.31b, an improvement in the resolution of the upper regions of the 
frequency spectrum is evident. In Figure 4.31b, a plot of the two-channel, four-subchannel 
decomposition of highpass test sequence (4.74b), the sequence component at digital frequency 
2-m-113/256 has been localized as distinct from scale zero. Furthermore, as compared with the 
two-channel, three-subchannel plot of figure 4.30b, the test sequence component at digital 


frequency 2:%:27/256 exhibits improved localization. 


TABLE 4.2--PARTITION OF SPECTRUM RESULTING FROM FIVE-STAGE, THREE-CHANNEL, 
ZERO-SUBCHANNEL MULTIRESOLUTION DECOMPOSITION. 


| Decomposition | Multiresolution Spectral bin Spectral bin 
lower bound —_ upper bound 


Psaieos [ys | os | 
Saal es | ca 


Scale 3.5 f/162 f./81 


Scale 4.0 f/243 f/162 


Scale 4.5 f/486 f./243 





As indicated previously, decomposition into three primary channels provides base three 


logarithmic partition of the frequency spectrum. Specifically, the first stage of decomposition 
divides the frequency spectrum into spectral bins covering the regions f,/6 to 2-f,/6 and 2:f,/6 to 


3-f£/6. During the second stage of decomposition, the lower third of the frequency spectrum Is 
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partitioned along intervals f/18 to 2:-f/18 and from 2-f/18 to 3-f/18. At stage m, the signal 
detail extracted will occupy the spectral regions f,/(2:3") to 2-f,/(2:3) and from 2-f,/(2:3") to 


3-f./(2:3"). Table 4.2 presents tabulated information regarding the spectral subdivision 


occurring from a three-channel, zero-subchannel multiresolution decomposition. 
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Figure 4.32a--Surface plot of decomposition of sequence (4.74a) using five-stage, three-channel, zero-subchannel 
multiresolution structure constructed from cascaded QMF banks. The decomposition was perfonned with 
pseudo-OMF filters designed in [11]. 


Finally, Figures 4.32a and 4.32b present the results of three-channel, zero-subchannel 
decomposition of the lowpass test sequence (4.74a) and highpass test sequence (4.74b), 
respectively. Because of the base three logarithmic partitioning of the frequency spectrum, 
localization of sequence spectral components in the lower portion of the frequency spectrum will 


be improved over the localization resulting from a two-channel decomposition. Specifically, for 
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sequence (4.74a), component at digital frequency 2:m:2/512 is located, according to Table 4.2, in 
the vicinity of the boundary between scales 4.0 and 4.5. In fact, most of the energy of this 
component appears in the channel corresponding to scale 4.0. Similarly, the component at 
2:%°9/512 is located in the vicinity of the boundary between scales 2.5 and 3. The energy of the 
appears mostly at scale 2.5. Finally, the energy of the component at digital frequency 
2:%°31/512 is indicated at scale 1.5. An unexplainable indication appears at a scale of 0.5 
between samples zero and 32 and also at approximately sample 96. From Figure 4.24, the filters 
from which the three-channel filter bank was constructed produced greater stopband error than 
the two-channel filter bank. Consequently, energy from the sequence appears to leak, producing 


less precise localization of components of the test sequence. 


With regard to Figure 4.32b which presents the three-channel, zero-subchannel of the 
highpass test sequence (4.74b), the structure appears to localize sequence components 
approximately as well as the two-channel, zero subchannel. Comparing the two-channel, 
zero-subchannel frequency partition diagram Figure 4.20 with the corresponding diagram Figure 
4.24 for the three-channel, zero-subchannel structure, the two methods provide very similar 


resolution in the upper regions of the frequency spectrum. Specifically, the three-channel 
structure divides the spectral region between f,/2 and f,/6 into two channels while the 
two-channel structure divides the region between f,/2 and f/4 into two channels. This similarity 
is reflected in the Figure 4.32b. The component at digital frequency 2:7%-27/256 appears as the 


largest range of peaks at a scale of two. The component at digital frequency 2-7%-113/256 is 
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indicated at a scale of zero along the edge of the surface. The indication of the component at 


2-%:55/256, however, is less well localized. Indications of this final component appear to be 


present at scales of 0.5 and 1.0. 


jana [auued) 
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Figure 4.32b--Surface plot of decomposition of highpass test sequence (4.74b) using five-stage, three-channel, 
zero-subchannel multiresolution structure constructed from cascaded QMF banks. The decomposition was 
performed with pseudo-QMF filters designed in [11]. 


Finally, the three-channel structure demonstrated in Figure 4.32a presents a noteworthy 
phenomenon. The surface in the higher scale (lower frequency) regions possesses a jagged 
texture. When the surface resulting from expansion by a filter is characterized by rapid 
variations such as are shown in Figure 4.32a, the filter is said to lack regularity. The appearance 
of this jagged quality is noteworthy because it occurs in the high-scale channels containing the 


slowly varying components. For the the two-channel structures, the surfaces in the regions of 
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higher scales presented relatively smooth appearances. The non-regular surface observed in 
Figure 4.32a increases the difficult in reading of the surface plot. The issue of regularity will be 


addressed again in Chapter V. 


To provide an overall assessment of the three-channel structure, it affords greater 
resolution at higher scales (lower frequencies) without incurring the increase in structure 
complexity involved in performing decomposition into subchannels. Compared to the 
two-channel, zero-subchannel structure, the structure provided somewhat improved resolution. 
However, the performance in the high-frequency regions Is not as great. Consequently, included 
among the issues in deciding whether to employ three channels or two-channels with 


subchannels is whether resolution at higher frequencies is critical. 
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V. BASIS FUNCTIONS AND FILTERS FOR MULTIRESOLUTION 
DECOMPOSITION 


A. INTRODUCTION 


Chapter IV developed the theory of multiresolution decomposition and demonstrated its 
equivalence to the cascading of QMF bank structures described in Chapter III. In Chapter IV, 
the wavelets and scaling functions used for multiresolution signal decomposition were addressed 
only the most abstract sense. In the present chapter the nature of these functions will be 
addressed in greater detail. Section B outlines some of the basic concepts of the theory of 
two-scale difference equations, the class of equations from which wavelets and scaling functions 
were, In Chapter IV, shown to be derived. In Section C, some of the properties of Daubechies’ 
orthonormal wavelets and scaling functions will be considered. Finally, Section D will return to 
a filter bank perspective for the study of multiresolution decompositions. Some basics of the 


theory of perfect reconstruction and pseudo-QMF bank filters will be addressed. 


It must be noted that a number of approaches exist for the construction of wavelets and 
scaling functions which will not be addressed in the present study. One of the earliest methods 
which has subsequently become relatively well-developed ts to construct wavelet and scaling 
function filters from polynomial splines [36, 37]. In the construction of wavelets from splines, a 
polynomial spline is used for the scaling function. Since polynomial splines of non-zero order 
are not orthogonal with respect to integer translation, the scaling function family constitutes a 
frame in the sense described in Chapter 2. The biorthogonal basis functions which complement 


the polynomial splines to complete the frame can be constructed from, among other methods, 


using recursive (IIR) filtering techniques. Finally, the design of wavelets has been approached 


from an optimization theory perspective, as well [38]. 


B. TWO-SCALE DIFFERENCE EQUATIONS 


In equation (4.23) (repeated here for convenience), it was shown that scaling functions 
represent the solutions to two-scale difference equations: 
5° OCF) =Zhlk)- O(u-k). (5.19 

Before attempting to solve (5.1), it is first useful to consider some of its properties. The 
development which follows makes the assumption that 9 1s a real function. However, the results 
can easily be generalized to complex functions as well. First, integrating each side of (5.1) over 
Its region of support produces [27] 

=- fo(=) du=Z h(k)- J (uk) du. (52) 
Applying the changes of variables v = u/2 and w =u -k to the appropriate sides of (5.2) results in 


{ o(v) dv = E h(k) J o(w) dw. (5.3) 


cancelled on both sides of the equation resulting tn 


PAs l, (3.4) 


Two properties have now been specified for the families of scaling functions which 


represent solutions to two-scale equations. The first is (5.4) which indicates that the filter whose 
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impulse response is h(k) possess a DC gain of unity. The second, which was necessary to obtain 
(5.4) requires that 


imo Genacin 0. (5. 


Cs 


) 


In fact, if the objective is to construct orthonormal basis functions, it is desirable to specify that 


If the condition is imposed that the family { o(u - gy ae 7 constitutes an orthonormal 
basis, then another interesting property for h(k) can be obtained. First, the change of variables t 
= 2-u1s made for (5.1) resulting in 
G(t) = 2. Zh(k)- G(2-t—k). (on5) 


Now, tf 0 1s orthonormal with respect to integer translation, then j O(t -k) O(t -j) dt = }; . where 


> is Kronecker's delta function. Furthermore, J O(2-t -k) O(2:1 - J) dt = 6, 2 also holds. 
Consequently, evaluation of the projection of integer translates of 0 onto each other yields 


J O(t-—m)- o(t—n) di=4-2% h(k) - hj) - | o(2-(t—m) —j)-o(2- (t-—n)—k) dt. (5.7) 
J 


Applying the change of variables u = 2-(t - m) to (5.7) produces 


| o(t —m)- o(t — n) dt=2- LE h(k) -hG) J o(2-u—))-o(u+2-(m—n)—k) dt. (5.8) 
j 


Next, applying the consequences of the assumption of orthonormality to both sides of (5.8), 
results in 


Bina = 2° DEAK) HG) 8.42000. (5.9) 
J 


Finally, by the summation sifting property of the Kronecker delta function, 


i! 


J 


Equation (5.10) provides the basis for at least two conclusions. First, for the case when 


E hk) =>. (nin) 


In other words, if the two-scale difference equation solution 0 is orthonormal with respect to 
integer translation, then the sum of the squares of the coefficients of the corresponding 


generating filter is one half. 


A second interpretation, made by Vetterli [25], state that the impulse response h(k) of the 
generating filter is orthogonal respect to even translates of itself. This is equivalent to another 
interpretation. Define the sequence g(n) from the convolution of h(n) with its own time 
reversal h(n): 

g(n) = 2+ E h(k)- h(n —k). (5.12) 
If g(n) is defined in this manner, then (5.10) 1s equivalent to the expression 

gain) = Ein) iS oa (5.13) 

Now, (5.13) 1s equivalent to (4.78). Therefore, the filter g(n) 1s an a trous filter in the sense 
described by Shensa [31] or a half-band filter in the sense described by Mintzer [32]. 
Furthermore, (5.13) requires that the length of h(k) be an even integer [25]. If the length of h 
were an odd integer L, then g(L - 1) = 2:h(0)-h(L-1) #0, and (5.13) would not hold. Finally, the 
evaluating the Z-transform of (5.12) results in 


G(z) = H(z) - H(z) = H(z): H(1/z). (5.14) 


2 


By (5.14), H(z) is simply a spectral factor of a half-band filter. This concept will be employed 


later in Section D of this chapter within the context of designing filters for QMF banks. 


Solving for scaling equations given an appropriate filter can be performed in either the 
time domain or in the Fourier domain. The Fourier-domain solution begins with Fourier 


transformation of (5.1), which was shown in (4.36b) to produce an expression equivalent to 
9(@) = H(ei?) - (2). (5.15) 
By induction, (5.15) produces [27] 


$(@) = (0) - TL H(e!®”*). (5.16a) 


Now, if the family of scaling functions 6 is orthonormal, then J o(u) du = | which is equivalent 


to 6(0) —simeineterore, (>.16a) becomes 

$(@) = TH(e mae) (5.16b) 
The Fourier transform of the scaling function 0 1s simply the infinite product of progressively 
more dilated versions of the Discrete Fourier Transform (DFT) of the generating filter h(k). 
Moreover, as the product index k tends towards infinity, the DFT term H(e! a2") on the 
right-hand side of (5.16b) becomes infinitely spread out in the Fourier domain. The Fourier 
transform of >, therefore, has an infinite region of support which indicates the possibility of 
finite support in the time domain. Compact support of the scaling function @ 1s highly desirable 


since, because of how h(k) was defined in equation (4.23), only a compactly supported scaling 


function will produce a finite impulse response filter. 
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To solve the scaling function in the time domain, it Is first convenient to absorb the 
factor of two on the right hand side of (5.6) into the filter coefficients by defining hy,(k) = 


2-h(k). Consequently, an equivalent form for (5.6) becomes 


O(t) = 2 hyay(k) Oeste (Salas 


The time-domain solution of (5.17) begins by solving for values of the function 9 at integer 


points in its domain [39]. Consequently, the equation to be solved becomes 


p(n) = {he (kK) G(s a (5.176) 


Now, if the scaling function @ is strictly supported on the region [0, 2-L - 1] and its generating 
filter h(k) is of length 2-L - 1 for some inter L, then (5.17b) can be expressed in matrix form as 


hey(0) 0 0 = 0 


(0) (0) 
hey(2) hey (1 hey (0 
(1) a : " ) oN ) | | (1) 
; 0 «su hey (20) aoe 2) ee 2 
Mi 2ecbuael) 0 - 0 0 ha(2" bee Dey ea 
(5.18) 


It is obvious that (5.18) is simply an eigenvalue problem. The vector containing the values of 
the scaling function o for integer arguments is simply the eigenvector of a unit eigenvalue for the 


double-nght-shift matrix of filter coefficients in (5.18). 


Now that the value of the scaling function @ has been solved for integer arguments, the 


next step is to rewrite (5.17b) in a form that will yield half-integer values. 


O(5) = Zhay(k)- (nk). (5.18) 
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Equation (5.18) is obviously a simple convolution 


expressing it in matrix form: 


(0) 0 

o(1) 0 

o(5) (2) 9(0) 
(>) > ! ‘ 
o(2-L-1) o(2-L=3) © 
— g( 2h) OHO 2) eos 

0 0 


, nevertheless, insight is gained from 


0 

0 

: hy2)(0) 

hy2)(0) 

5.19 
(0) | o 
O(1) hy) (2eele= | ) 

a (2: 5 


From (5.19), a pattern begins to emerge. Extended (5.18) to the general dyadic case produces 


O( seer) = Zhay(k) (ge — Kk). 


Finally, adopting the vector notation 


ale 


m 


gm =( 


th 


oe 


b=( he) halt) 





oa 
— 


(5.20) 


2™.(2-L=1) 
2m 


: 


T 
ha(2:-L=1) 


and 
(<=) 0 
o(s4) 0 
o22) 
pin | Cae) GR) 
o( FE) oS) o(2) 
0 g( SER) ok) 
0 0 (FS) 


> 


equation (5.20) can be expressed in a more compact form: 


ge =O h. (S.21a) 


Similarly, if the vector y'”’ for the wavelet function y(t) is defined in a manner analogous to 


'""” for the scaling function $(t) and the wavelet generating filter g in in a manner analogous to 


the scaling function generating filter h, then, as a consequence of (4.51c), the values of w(t) can 


be obtained at dyadic points from 


yr =—™ g. (5.21b) 


Sample value s(n) 


Sample sequence number, n 


Figure 5./a--Function generated by four-stage dyadic expansion of a 20-point lowpass filter with cutoff frequency 
of 0.375-f,. The function represents an example of a function which lacks regulanty. 


The recursive operations described by (5.21a) and (5.21b) represent dvadic expansions of 


the functions o(t) and y(t). At each iteration the discrete sequences characterized by the vectors 


o'” and y'™ are expanded by a factor of two and then convolved with filters h(n) and g(n), 


respectively, whose coefficients are contained inh and g. After m iterations, the values of the 
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of the scaling function o(t) and wavelet y(t) will be known at all points t, in region of support 


such that t, = leo. 


A final issue which must be addressed within the scope of two-scale difference equations 
1s the topic of regularity. Regularity is defined to be the degree to which a recursively expanded 
function, such as a two-scale difference equation, converges to a continuous smooth function 
[26]. The conditions under which the regularity of a function can be guaranteed of themselves 
constitute grounds for an extensive study. However, one apparent example of a requirement for 
regularity was observed during experimentation performed during the course of the present 
study. For regularity, it is necessary that, for the filter to be expanded, the DFT exist primarily 
within the lower half of the frequency spectrum below the Nyquist frequency. Figures 5.la and 
5.1b show two examples of the dyadic expansion of two candidate filters for the generation of a 
scaling function. Both filters were designed by performing spectral factorization of lowpass 


filters designed using the window method [12]. For the first filter the cutoff frequency was 


selected as 0.375-f, and the second filter was assigned a cutoff frequency of 0.25-f,. As can be 
seen in Figure 5.1a, the function possesses an erratic appearance. As the bandwidth of the 
candidate function generating filter increases, the generated function becomes less smooth and 


more and more erratic. 


On the other hand, Figure 5.1b shows a function generated through dyadic expansion of a 


filter designed with exactly the same specifications except that the cutoff frequency was 


specified as 0.25-f,. The plot of the function generated by the second filter displays significantly 
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greater smoothness. Figure 5.1b represents an example of a function which, as far as can be 


determined from the plot, possesses the quality of regularity. 


sample value s(n) 





Sample sequence number, n 


Figure 5./b--Function generated by four-stage dyadic expansion of a 20-point lowpass filter with cutoff frequency 
of 0.5-f,. The function represents an example of a function which possesses regularity. 
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C. DAUBECHIES' FAMILY OF ORTHONORMAL WAVELETS 


Ingrid Daubechies of Bell Laboratories madé a landmark contribution to the field of 
multiresolution analysis when, in 1988, she published a mathematical paper on orthonormal 
wavelets [21]. Aside from providing a general overview of the progress which during recent 
years had occurred within the field of wavelets, Daubechies' paper described the construction of 
a class of functions which possessed many of the most desirable aspects of wavelet functions. 
The functions were specified in terms of their generating filters. Even-length filters from length 
four to 20 were tabulated in the paper. Each of the filters satisfied (5.4) and (5.11). 
Furthermore, they each maximized the number of zeros on the unit circle at z = -1, a condition 


which ensures maximum regularity [26]. In this section plots of two wavelet-scaling function 


pairs will be presented, and one will be examined more thoroughly from the perspective of 


implementation of the multiresolution structures from Chapter IV. 
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Figure 5.2a--Superimposed plots of Daubechies' orthonormal wavelet and scaling function on (0, 3]. Plots were 
generated by four-stage dyadic expansion of corresponding filter listed in (21). 





Figure 5.2b--Polar plot of zero locations for generating filter for Daubechies' orthonormal wavelet and scaling 
function supported on (0, 3]. 
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Figure 5.3--Two-channel, two-subchannel multiresolution decomposition of 256-point lowpass test sequence 
(4.74a). Decomposition was performed using generating filter for Daubechies' orthonormal basis function 
supported on (0, 3]. 


Perhaps the most famous of ideal orthonormal scaling function-wavelet pair is 
plotted in Figure 5.2a. The functions are generated from a four-point filter and are supported on 
(0, 3]. As can be inferred from viewing the plots, the smoothness of these functions 1s marginal. 
Figure 5.2b provides a polar plot of the zeros of the generating filter. Of its three poles, two fall 
on the unit circle. Finally, as indicated in Chapter IV, the wavelet transform, which essentially 
consists of the detail function (4.33), is comprised of a plot of weighted, shifted versions of the 
wavelet. Consequently, it is expected that the texture of a surface plot of a wavelet 
decomposition should reflect the shape of the wavelet from which it is constructed. Figure 5.3 
supports this expectation. In particular, the peaks in the which appear at higher scales reflect the 


shape of the wavelet function plotted in Figure 5.2. 
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Figure 5.4--Superimposed plots of Daubechies' orthonomnal wavelet and scaling function on (0, 13]. Plots were 
generated by four-stage dyadic expansion of corresponding filter listed in [21]. 





Figure 5.5--Polar plot of zero locations for generating filter for Daubechies' orthonormal wavelet and scaling 
function supported on (0, 13]. 


Figure 5.4 presents superimposed plots of Daubechies' orthonormal basis function and 
scaling function supported on [0, 13]. This pair was used extensively throughout Section IV.E 


to illustrate the performance of multiresolution decomposition techniques. The zeros of the 


14] 


generating filter are plotted in Figure 5.5. Of the 13 zeros of the generating filter's 


characteristic polynomial, seven occur at the location z = -1. 
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Figure 5.6--Time-domain demonstration of perfonnance of QMF bank constructed from generating filter for 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Left-hand axis indicates equivalent impulse 


response while right-hand axis indicates deviation t(n) - 5(n - N,,,,) from ideal perfect reconstruction 
performance. 


In Section IV.E, multiresolution decomposition techniques were addressed from the filter 
bank perspective. In Chapter III, the degree to which specific examples of filter banks satisfied 
the perfect reconstruction criterion and the power complementary property was demonstrated. 
The question, therefore, arises as to the degree to which the generating filters for Daubechies' 
orthonormal wavelets and scaling functions satisfy the perfect reconstruction and power 
complementary properties. Figures 5.6 and 5.7 illustrate the performance of these functions 
employing the same plots used in Chapter III. Figure 5.6 presents a plot analogous to Figure 
3.14. Upon comparing Figure 5.7 with Figure 3.14, it is apparent that Daubechies’ generating 


function provides vastly superior performance for filter bank applications. The filter bank 


demonstrated in Figure 3.14 deviated from the ideal response of a perfect delay by 2*10 " 
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Daubechies' 14-point generating filter deviates from an perfect delay by 1.4«x10 This 
represents six orders of magnitude in performance improvement over the time-domain 


performance of the QMF bank considered in Section IILC. 


Figure 5.7 indicates similar results for the frequency-domain performance of a QMF 
bank constructed from Daubechies' 14-point generating filter. The desirability ofa 
power-complementary pair of analysis and synthesis filters was addressed in Chapter III. The 
QMF bank of Figure 3.15 deviates from being a perfectly power complementary pair by 10° dB. 


A QMF bank constructed from Daubechies' 14-point generating function deviates from being a 


perfect power complementary pair by 1.2x 10°" dB. This represents ei ght orders of magnitude in 
improvement. Consequently, Daubechies’ generating function for wavelets and scaling functions 
supported on [0, 14] provide for outstanding QMF bank performance conforming very closely to 


the power complimentary and perfect reconstruction properties. 
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Figure 5.7--Demonstration of the degree to which generating filter for Daubechies' orthonormal wavelet and 
scaling function on [0, 13] produces a power complementary pair for QMF bank purposes. The left-hand axis 
indicates the frequency responses for the filter and its modulated complement. The nght-hand axis indicates the 
degree to which the filter pair forms a power complementary set. 
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D. DESIGN OF FILTERS FOR OMF BANKS 


[In the final Section of this chapter, procedures for the desizn of filters for filter bank 
applications will be examined. The steps for the design ofa filter set for a two-channel QMF 
bank will be outlined. After that, the process of a filter design for a four-channel pseudo-QMF 
bank--a filter bank which only approximately satisfies the perfect reconstruction and power 
complementary properties--will be outlined. Procedures for the design of perfect reconstruction 
banks of an arbitrary number of channels exist [1 1, 34], but will not be addressed in the present 


study. 


To construct a two-channel filter bank with real-coefficient filters which possesses the 
perfect reconstruction property, it 1s necessary, as was shown by (3.20a), that the Z-transforms 
of the synthesis and analysis filters, both of order N, for a given channel k be related by [13] 

F.(z)=z-H, (1/z’). (5122) 
Equation (5.22) 1s equivalent to the notion that impulse response of the synthesis filter 1s equal 
to the time reversal of the impulse response for the analysis filter. Additionally, from (3.19), the 
equivalent transfer response T(z) of the filter bank structure must, to satisfy the perfect 


reconstruction criterion, be of the form 


T(z) =4-[Fo(z)-H(z) + F,(z)-H,(z)] =z’. (5.23) 
For convenience, an equivalent transfer response G,(z) for the branch of the filter bank 
corresponding to each channel can be defined: 
Gi(Z) = 5° Hx (Z) > F(Z). (5.24) 
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Now, assume that G,(z) represents the transfer function for a half-band filter. If G,(z) 


satisfies the half-band filter property in the sense described by Mintzer [32], then the impulse 
response g,(n) corresponding to G,(z) satisfies [13] 

(g.)1 (1) = Sy (5.25) 
where N represents the orders of the filters h,(n) and f,(n). Therefore, by (5.25), the only 
non-zero, even-numbered element of the impulse response of g,(n) 1s the sample g,(N). 
Continuing, the objective is to satisfy (5.23), or equivalently, 

T(z) = G,(z) + G,(z) =z". (5.26) 
One way to satisfy (5.26) is to specify that the deat bee elements of g,(n) each be of the 
opposite sign, or 
eee 1) = 22 © 1). (27a) 
Because of (5.25), (5.27a) 1s equivalent to 


g,(n) = (-1)"g,(n). (3270) 
But, (5.27b) is equivalent to 


G,(z) = Gol-z). Gr27c) 


From the preceding development, it may be inferred that H,(z) and F,(z) are 
equal-ordered spectral factors of a halfband filter G,(z). Furthermore, g,(n) is simply a 
modulated version of g,(n). Therefore, the frequency response Ge ”) consists simply of the 


frequency response of G,(e ”) shifted in the frequency domain by x. Finally, since G,(z) 1s a 
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polynomial with two equal-ordered spectral factors which are simply the reversals of each other, 
certain constraints on the locations of the zeros of G,(z) arise. First, for each zero z, of G,(z) on 
the real axis corresponding to H,(z), there must be a reciprocal zero l/Z, on the real axis 
corresponding to F,(z). Secondly, since it has been specified that the filter impulse responses 
h,(n) and f,(n) be strictly real, the zeros of G,(z) must occur in complex-conjugate pairs. 
Furthermore, for every zero Z4 located inside the unit circle corresponding to H,(z), there must 


be a corresponding reciprocal zero 1/z, corresponding to F,(z) located outside the unit circle. 


Consequently, all zeros not on the unit circle or the real axis must occur in complex conjugate 
quadruples. Finally, in order for H,(z) and Bue) to be of equal order, for each 
complex-conjugate pair of zeros z, and 7 located on the unit circle corresponding to H,(z), an 
identical set corresponding to F,(z) must occur in the exact same location. In other words, the 


zeros on the unit circle for G,(z) must occur as double zeros. Consequently, G,(z) must be a 


filter with a strictly real, non-zero, symmetric, zero-phase filter. 


The requirements for the locations of the zeros for G,(z) suggest a design technique for 


filters for perfect reconstruction filter banks. The first step 1s to design a zero-phase, symmetric 


half-band, lowpass filter G,(z). This may be accomplished by any standard method such as the 


McClellan-Parks routine or by the window method [12]. Such a filter will be symmetric and 
will satisfy the half-band filter criterion (5.25). In order for G,(z) to be strictly non-zero, the 


maximum stopband ripple €, must next be measured. If this mpple €, is added to the central 
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element of the filter impulse response g,(n), the result will be an impulse response of a modified 
filter 
Sog(N) = ¥o(n) + Ep'd(n - N) (5.28) 
with a corresponding transfer function 
Cee) — Soe) Ep. (5.28b) 


Because of the addition of the term €, the frequency response of G,,(z) will be strictly non-zero. 


h (n) 


eee 
“0.3 





Discrete time, n 


Figure 5.8--Impulse response of 30-point, lowpass filter designed for QMF bank. 


The next step requires the rooting of Gp,(z). The zeros of Gop(z) must be grouped into 


three subsets: zeros located inside of the unit circle, zeros located on the unit circle, and zeros 
located outside of the unit circle. The zeros inside of the unit circle belong strictly to H,(z) and 


the zeros outside of the unit circle belong strictly to F,(z). The zeros on the unit circle should all 


be double zeros occurring in sets of complex conjugate pairs. One set of each double complex 


conjugate pairs belongs to H,(z) while the other belongs to F,(z). If the zeros corresponding to 
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H(z) are re-expanded into a polynomial, the coefficients of the result should be strictly real and 
belong to h,(n). 
As an example, a 30-point QMF filter was designed whose impulse response is plotted in 


Figure 5.8. The halfband filter was generated from applying a Kaiser window to a 59-point sinc 


function 


Sin) —-(n=N) 
go(n) = sli) ) 


After adding the peak stopband error €, to impulse response element g,(30), the corresponding 
polynomial was rooted and zeros grouped as was described above, resulting in a polynomial 


whose coefficients represent the impulse response ho(n) of a QMF-bank lowpass analysis filter. 
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Figure 5.9--Time-domain demonstration of performance of QMF bank constructed from the 30-point lowpass 
filter designed by spectral factorization of a half-band filter produced by the window method. Left-hand axis 


indicates equivalent impulse response while right-hand axis indicates deviation t(n) - 6(n - Ny,,) from ideal 
perfect reconstruction performance. 
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Figure 5.10--Demonstration of the degree to which 30-point lowpass filter designed by spectral factorization of a 

half-band filter produced by the window method produces a power complementary pair for QMF bank purposes. 

The left-hand axis indicates the frequency responses for the filter and its modulated complement. The right-hand 
axis indicates the degree to which the filter pair deviates from a power complementary set. 


Figures 5.9 and 5.10 indicate the results of application to the filter the analysis methods 


demonstrated in Chapter III for filter banks. Figure 5.9 plots the equivalent transfer response of 
a filter bank constructed from h,(n), and its deviation from the ideal response of an ideal delay. 
As can be seen from Figure 5.9, this method works quite effectively. The peak deviation from 


an ideal delay was less than 3x 107°. This result compares favorably with the response for the 
generating filter for Daubechies' wavelets examined earlier. The power complementary 
properties are demonstrated in Figure 5.10. This filter did not quite perform as well as the 


Daubechies'’ wavelet generating filter performed. The peak deviation from a response for a 


“ll for the 


power complementary pair was approximately 6x 10° compared with 1.2*10 
Daubechies filter. This result is also slightly inferior to the performance of Daubechies' wavelet 


generating function by roughly two orders of magnitude, but still significantly superior to the 


QMF bank filters considered in Chapter III. 
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Finally, within the context of a study addressing wavelet transforms the question arises as 


to the use of perfect reconstruction QMF filters for generating wavelets. First of all, the filter 


impulse response h,(n) does not strictly satisfy (5.11). Specifically, 


~ h2(k) = ++7.5x 10%. 
: 2 


Consequently, any scaling functions or wavelets generated by h,(n) would not be orthonormal. 
Additionally, if a two-scale difference equation for h(n) is set up in the form of (5.18), the 


nearest eigenvalue to unity for the double-right-shift matrix of coefficients of h,(n) is 1.0052. 


As aresult, only an approximate system of a two-scale difference equation can be constructed 


from h(n). Based on the preceding two measures it may be concluded that although, as was 
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Figure 5.11--Plots of "pseudo-scaling function" and "pseudo-wavelet” generated from four-stage dyadic expansion 
of 30-point filter designed from spectral factorization of a half-band filter. 


shown in Chapter IV, wavelet decompositions are exactly equivalent to structures of cascaded, 
perfect reconstruction QMF banks, the converse is not true. In fact, in the present example, the 


perfect reconstruction filter bank only approximately satisfies the properties described for 
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orthonormal wavelet decompositions. Nevertheless, if these issues are ignored and a dyadic 
expansion of h,(n) is performed, the resultant functions, after four iterations, appear as plotted in 
5.11. The functions possess obvious smoothness and present appearances similar to those of 


Daubechies’ wavelets. 


To generate a filter bank for a four-subchannel multiresolution decomposition, a 
pseudo-QMF bank was designed using a spectral factorization technique similar to that used for 


the preceding example. The strategy was to approximate the power complementary property. 
Using the McClellan-Parks algorithm, two filters h,(n) and h,(n) were designed whose 
passbands covered the regions [0, 1/4] and [1/4, 1/2], respectively. The upper half of the 


frequency spectrum below the Nyquist frequency was covered by modulating the first two 


filters: 


h,(n) = (-1)"h,(n) 


and 

h,(n) = (-1)""ho(n). 
The design of the first two filters involved an iterative process in which, at each stage, the 
closeness of the filter set to the power complementary process was checked. At each step, the 
sum of the frequency responses was adjusted by modifying the boundaries to the passband and 
stopband for each filter. Tables 5.1a-c document the iterations. Each column of Tables 5.1a and 
5.1b indicates the frequency at which the desired gain in the top row Is to be met for the filter to 
be factorized spectrally. Table 5.lc indicates, at each iteration, the peak deviation in the 


transition bands between each filter. Based on the deviation in the transition bands between 
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ho(n) and h,(n) and between h,(n) and h,(n), the filter indicated by h)(n) was not modified after 


the third iteration. Eight iterations were necessary to obtain acceptable results for the transition 


band between h,(n) and h,(n). 


TABLE 5.1A--RECORD OF ITERATIVE ADJUSTMENTS TO PASSBAND AND STOPBAND FREQUENCY 
BOUNDARIES IN DESIGN OF FIRST FILTER H,(N) FOR FOUR-CHANNEL, PSEUDO-QMF BANK. 
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TABLE 5.1B-RECORD OF ITERATIVE ADJUSTMENTS TO PASSBAND AND STOPBAND FREQUENCY 

BOUNDARIES IN DESIGN OF SECOND FILTER H,(N) FOR FOUR-CHANNEL, PSEUDO-QMF BANK. 
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Upon review of Table 5.1c and of Figures 5.13 and 5.14 it becomes obvious that the 
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Iteration 4 







Iteration 5 






Iteration 6 






Iteration 7 






Iteration 8 





performance obtained from the four-channel filter bank is nowhere near what was obtained by 


any filterbank considered at any other time during this study. The peak delay of the equivalent 


impulse response for the filter bank was of the order of magnitude 5x 10°. Furthermore, the 


eZ 


peak deviation from the power complementary standard could only be reduced to tenths of 


decibels. 


TABLE 5.1C--TRANSITION BAND DEVIATION FROM POWER COMPLEMENTARY, BY ITERATION, 
FOR DESIGN OF FOUR-CHANNEL PSEUDO-QMF FILTER BANK. 





Peak deviation in transition 
bands between h,(n) and h,(n) Peak deviation tn transition 
and between h,(n) andh,(n) band between h,(n) and h,(n) 


Iteration 1 Excessive Excessive 


teration3[ 007—S«dYSCS 





Discrete time, n 


Figure 5.12--Impulse responses hp(n) and h,(n) of filters designed by spectral factorization of fourth-band filters 
designed using McClellan-Parks algorithm. 
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Figure 5./3--Time-domain demonstration of performance of four-channel QMF bank constructed from filters 
designed from spectral factorization of fourth-band McClellan-Parks filters. Left-hand axis indicates equivalent 
impulse response while nght-hand axis indicates deviation t(n) - 6(n - N,,,,) from ideal perfect reconstruction 
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Figure 5.14--Demonstration of the degree to which four-channel, pseudo-QMF bank designed by spectral 
factonzation of a McClellan-Parks fourth-band filters produces a power complementary pair for QMF bank 
purposes. The left-hand axis indicates the frequency responses for the filter and its modulated complement. The 
nght-hand axis indicates the degree to which the filter pair deviates from a power complementary set. 
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VI. EVALUATION OF MULTIRESOLUTION TECHNIQUES FOR 
DETECTION APPLICATIONS 


A. INTRODUCTION 


In this chapter, additional performance characteristics of the multiresolution structures 
developed in Chapter IV will be examined within the context of evaluation of suitability for 
detection applications. In Section B, the computational efficiency of multiresolution structures 
will be considered. In the Section C, receiver operating characteristics will be plotted for 
various multiresolution structures applied to the test sequences generated by (4.74b) and (4.74a). 
Finally, Section D will demonstrate resolution of scale within the context of steady-state 


harmonics and chirped signals. 


The basis for evaluation of multiresolution performance will be the spectrogram, one of 
the most commonly used techniques for detection and spectral estimation. Figures 6.1a and 
6.1b present the spectrogram decomposition of the two 256-point test sequences generated by 
(4.74a) and (4.74b), respectively. The plots were generated by plotting the frequency spectrum 
below the Nyquist frequency for a sequence windowed with a 128-point Hanning window 
shifted in increments of two samples. This representation was selected because of its frequency 
resolution and because it presents a similar number of points compared to the multiresolution 


decompositions presented in Chapter IV. 


In Figure 6.1a, the lack of frequency resolution in the extreme lower range of the 


frequency spectrum is apparent. For lowpass test sequence (4.74a), the spectral components at 


oD 


digital frequencies 2:%-2/512 and 2:2:9/512 are not clearly resolved. The component at 
2-m:2/512, which appears at spectral bin location 0.5, because of its relatively greater power, 
spreads into adjacent spectral bins. Consequently the combenent at 2:7-9/512, which coincides 
with spectrogram bin 2.25, is partially masked by leakage from the component at 2:m:2/512. The 


only indication of the component at 2:%-9/512 is evident on the leading edge of the surface plot. 


7 < >. 
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Figure 6.1a--Spectogram representation of lowpass test sequence (4.74a). The plot was generated using a 
|28-point Hanning window shifted in increments of two samples. Only the half of the frequency spectrum below 
the Nyquist frequency ts plotted. 


On the other hand, Figure 6.1b provides an obvious indication of the advantages of 
constant spectral bin density: The three harmonics of (4.74a) are quite easily resolved. The 


plots have the additional advantage of a smoother presentation. Spectrogram decompositions 
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smooth out the phase information while multiresolution presentations indicate individual 


oscillation cycles of the sequences decomposed. 
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Figure 6.1b--Spectogram representation of highpass test sequence (4.74b). The plot was generated using a 
128-point Hanning window shifted in increments of two samples. Only the half of the frequency spectnim below 
the Nyquist frequency is plotted. 


B. COMPUTATIONAL EFFICIENCY OF SIGNAL ANALYSIS TECHNIQUES 
Computational efficiency comprises one of the principal advantages ascribed to Fast 

Fourier Transform (FFT) methods of digital signal processing. The required number of 

multiplications provides a common number for evaluating digital processing techniques. In the 


case of FFT decomposition of a sequence of length M, the required number of complex 


multiplications N, is expressed by 


Nx == log,(M). (6.1) 


ro 


Furthermore, application of the window sequence to each segment of a sequence of length L 


requires M multiplications. Finally, if the window is shifted in increments of k, then the FFT 


will be applied to =“+1 sequences. Therefore, the total number of multiplications necessary 


to generate a windowed FFT decomposition 1s 
Nx = (se +1 Jog, (M). (6.2a) 


In general, the multiplications performed in FFT decomposition represent complex 
multiplications. Evaluation of a complex multiplication requires four real multiplications and 
two additions. Consequently, assuming that real data and window sequences are used, the 


number of multiplications for the FFT decompasition is actually 

Nx = 2( E+ 1) Me log,(M). (6:2b) 
For each of the decompositions plotted in Figures 6.la and 6.1b, the sequence length L = 256, 
and the window length M = 128. Furthermore, the window was shifted in increments of k = 2 
samples. Consequently, using (6.2b), the number of real multiplications N.. required to perform 
the decomposition was N,, = 14,909,0440. 

For the multiresolution decompositions of Chapter IV, Section E, FF T-type radix-2 

algorithms are not available. Consequently, each stage of decomposition requires the 
convolution of two real sequences. Therefore, for decomposition of sequences of length L,, 


using filters of length M, for decomposition into two primary channels, the required number of 
real multiplications for the first stage is 


N, = 2:M,'Lo. 
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After convolution by the analysis filters of each stage, the approximation and detail sequences 


are both decimated, providing sequences of length 


”" Lot+Me-l 
= a. 


ra 


Ly 


By induction, the length L, of the approximation and detail sequences at stage p is given by the 


geometric series 


Lo ea] 
Lop = ge t(Me- 1): 2 ag 
ce q= Cd 


= BL eee | (6.3) 


2P 


= $e (Lo-(Me-1)) +(Me=1) 


Since, at each stage, each of two sequences is convolved with a filter of length M,, the total 


number of multiplications N, necessary to perform a K-stage multiresolution decomposition into 


primary channels 
K-1 
Neacele- 2 ba. (6.4) 
p=0 


Substitution of (6.3) into (6.4) produces 


K= 
Ne =2-K+(Me=1)+2- Me: (Lo-(Me-1) ); E sr. (6.5) 
a 
Finally, evaluating the geometric series in (6.5) produces 


Equation (6.6) covers all of the multiplications necessary to perform the decomposition into 


primary channels in accordance with block diagram 4.19. 


lo 


If the primary detail channels are each to be decomposed into J subchannels using a 
subchannel filter bank constructed from filter of length M,, additional computation is necessary. 


First, if a sequence of length L, 1s convolved with a filter of length M, and subsequently 


-°? of the resultant sequence will be 


decimated by a factor of J, the length L 
(p) 
cf = (Looe ee 
The number of multiplications N:P to perform this signal reduction at stage p is expressed by 
N,? = M,'L,, 

which by substitution of (6.3) becomes 

(p) 1 ; 

 =Me-(4-(Lo-(Me-1))+(Me-1)), 
If there are a total of J subchannels per stage and a total of K stages, then the total number of 


multiplications N_ required to perform all subchannel decomposition operations is evaluated as 


EN? 


p=1 


Ny 


JK +My (Me=1)+J-Mg-(Lo-(Me-1) ); 5 ee . (6.7) 
p=i° 


J+ K+ My: (Me= 1) +J- Ms: (Lo -(Me~1)}-(= 35) 


If the decomposed sequence of length io ’ is expanded by a factor of J and subsequently 
convolved with a filter of length M., the length Le of the resultant sequence 1s expressed by 
(P) — (p) = 
Legh =oJ-LePaevh, aoa vies): 
The number of multiplications N,'” necessary to perform this operation for J subchannels at 


decomposition stage p is therefore 
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N,? =J-M,L,” = J-M,-(L, + 2-(M, - 1)). 
Furthermore, the total number of multiplications necessary to expand all subchannel 


decompositions 1s expressed as 


K | ip) Ky | 
Ne =J-M,- 2 Le =J-Ms E | ae (Lo-(Me~1) )+(Me- 0) (6.8) 


eS 


Substituting (6.3) into (6.8) and evaluating the series produces 


Ne 


es 2 / “(M.-1)1 
J-M, 2 LP (Le (M, 1) )+(M. 2 (M, - 


JK My: ((Me 1) #2: (Ms-1) ) +5. M, (Lo-(Me-1))- E 5 . (6.9) 


K+ Ms: ((Me=1)+2:(Ms=1) )+J-My: [Lo -(Me =) (I=) 


To calculate the number of ti adene necessary to re-expand a detail component © 
sequence reconstructed from decomposition into subchannels, the notation one will be 
employed to indicate the length of the sequence at a stage of interest. The superscript ‘”’ 
indicates the number of the stage of decomposition at which the detail component sequence was 
extracted. The subscript , denotes the stage of reconstruction. The stages of reconstruction are 
labeled in decreasing order from p-1 to 0 for a component extracted at the a stage of 
decomposition. At stage p, the length of the detail sequence, after reconstruction from 
subchannel decomposition, 1s given by ee = oe Now, upon expanding the detail sequence of 
length ee ’ and convolving it with the channel synthesis filter of length M., at stage p-1, the 
expanded sequence which 1s obtained has length 


LP = 2," + M,-1=2-L,+4(M,-1)+M,- 1. (6.10) 
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Inductively, after q stages of decimation followed by convolution with the appropriate channel's 


synthesis filters, (6.10) becomes 


- | 
Low? = 29. [Lp #2-(Ms=1) }+ (Me 1) £2 


24. (Ly #2-(My=1) )#(Me=1)-(24=1) | (6.11) 


29. (L, +2 (M,—1)+(M.- 1) )=(Me= 
Now, substituting L, from (6.3) into (6.11) produces 


Lea = 2%: (a5: (Lo-(Me~1) )+(Me=1) +2-(Me= 1) +(Me-1) )-(Me=1) 
(63) 
Consequently, the total number of multiplications required for a full expansion by channel of a 
decomposed sequence can be evaluated by summing up the contributions from each convolution 
performed in each channel. It stands to reason, that the number of multiplications News , 
involved in expansion of J detail component sequences (one detail component sequence for each 
of J subchannels) of length L,”, at stage p, to obtain a sequence of length L,,,'?’,at stage p - 1, is 
indicated by 
N My = IM,L,.,”. 


a= 


‘P) necessary to completely 


Inductively, therefore, the total number of multiplications N, 
expand the J decomposed sequences in channel p is obtained from 


Pp 
N,” =] - Meow Logie (6.13) 
q=l 
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Therefore, to obtain Nae It is necessary to substitute (6.12) into (6.13) and evaluate the series, 


which produces 


p 
nor - J Me-| $e: (Lo (Me ~1) ) #2: (Me-1)42:(Ms—1)] & (29)—p-J- Me: (Me= 1 


q=! 


il 
bo 
Cony 
< 


2 [+ (Lo-(Me= 1) +2: (Me-1)42-(My=1)] QPP Me (MoD 


T 
= 


|  (Lo=3-(Me=1)=2:(My=1) )#2-J- Me: (2: (Me=1)42-(M,-1))- 2 
-2:J-Me: (Lo-(Me~1) )-=p-J-Me-(Me~1) 
(6.14) 


And lastly, to obtain No, the total number of multiplications to compute the full 


expansion by channel, it is necessary to sum the contributions N,° ’ from each channel. This 


sum appears as 


No = 2 J-K-Me-(Lo-3+ (Me=1)-2-(Ms~1) | 
#2:J-Me-(2:(Me=1) +2: (Ms-1)): y 2° 
: — (6.15) 
-2-J-Me-(Lo-(Me~1))- E 


pen 


K 
-J-Me(Me-1)- 2 p 
p= 
The first series in (6.15) are simply the geometric series evaluated previously in this 


development. The last series can be solved using Bernoulli polynomials [40] which produce 


p=+-(K +k), 
= 


P 


Substituting each of the evaluated series into (6.15) and collecting similar terms produces, for 


No; 
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a 2-K-J-Me: (Ly -3-(Me~1)-2-(M,~-1)) 
#4:3-Me:(2-(M. —1)+2-(M :-1)) Ce 
-2-J-Mg- (Lo -(Me=1)}-(1= 4) 


—4 J) Me: (Me~1)-(K?+K ] 


(6.16) 


Equation (6.16) applies whether or not subchannel decomposition is performed. If 


subchannel is not performed, then the detail sequence at each stage is decomposed into J = 1 
subchannels and the filter g.(n) used for subchannel decomposition possesses as its impulse 
response g,(n) = 6, , where 6, , is Kronecker's delta function. Consequently, because of how 


g,(n) is defined, M=I and, in (6.16), each of the terms in involving M, - | become zero. 


TABLE 6.1--COMPARISON OF NUMBER OF MULTIPLICATIONS NECESSARY FOR MULTIRESOLUTION 
AND SPECTROGRAM SIGNAL REPRESENTATIONS. 










Primary Primary 
Channel Channel Subchannel Subchannel 
Decomposition Re-expansion- Decomposition Re-expansion Total . 


hae 


— | 












Two-Channel, 
Zero-Subchannel 


Two-Channel, 
| Two-Subchannel 





Two-Channel, 
-Three-Subchannel 


Two-Channel, 
Four-Subchannel 





Spectrogram 





Summarizing, there are four contributors which must be evaluated to calculate the 
number of multiplications N, necessary for a full channel by expansion of a decomposed 


sequence. The basic decomposition, N., is provided by (6.6). The decomposition into 
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subchannels, N. is calculated using (6.7). To re-expand subchannel decompositions, the number 
of multiplications is indicated by (6.9). And finally, for the full expansion by channel, the 
number of calculations 1s provided by (6.15). For some of the multiresolution structures 


demonstrated in Section IV.E, the number of multiplications is tabulated in Table 6.1. 


In considering Table 6.1, it must be noted that the decompositions listed do not represent 
the optimum representation. In the case of the multiresolution decompositions, the essential 
information about the decomposed signal is contained in the primary channel and subchannel. 
The re-expansion by channel merely provides a manner of displaying the information which is 


more convenient than plotting coefficients in a non-uniform lattice of the type of Figure 4.8. 


In the case of the spectrogram, shifting a 128-point window by increments of two 
samples provides what is likely excessive overlap. This degree of overlap was used in an 
attempt to capture, to the greatest extent possible, the time-varying features of the test sequences 
on which decomposition was performed. If, for instance, the 128-point window was shifted in 


increments of 64 samples to provide for a 50% overlap, the number of multiplications N, would 


be N,=688,128. This overlap, however, would provide poorer indication of the non-stationary 
features of the sequences for which decomposition was performed. A 50% overlap with a 
window length of 128 used to decompose a sequence of length 256 would provide for three time 


shifts. 


With regard to the multiresolution decompositions, excluding the calculation used to 


re-expand the components of decomposition, a significant advantage in computational efficiency 
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is evident. Because of the factor-of-two decimation performed at each stage of decomposition, 
the sequence length at each stage 1s one half the length of the sequence of the previous stage. 
Consequently, the number of calculations necessary to perform the calculation is related to a 


geometric series In one half. 


Cc. PERFORMIANCE OF DECOMPOSITIONS IN THE PRESENCE OF NOISE 


SS ee ee ee 


Many signal processing techniques are employed with the intent of extracting a signal 
from noise in which it is embedded. The spectrogram, one of the most commonly used signal 
processing tools, represents an attempt to localize spectral components of a signal. If the noise 
in which the signal 1s imbedded can be characterized as "white," possessing a flat power spectral 
density, a signal which can be localized in a een sense will be detectable in greater noise 
intensity than a signal which 1s not spectrally localized. Consequently, when attempting to 
detect a signal by means of spectral localization, the effective signal-noise ratio can be reduced 
to the ratio of signal power to the power of only the noise contained within the spectral region of 


interest. 


To detect a time-varying process, the detection process becomes a two-dimensional 
problem, requiring localization in both frequency and time. Within the context of 
time-frequency localization, multiresolution techniques present the potential for an alternative 
method of detection. Many of the advantages of multiresolution transforms which lend those 
decompositions to signal decomposition also provide the basis for extraction of embedded 
signals from noise. Multiresolution transforms provide good time resolution at high frequencies 


and frequency resolution at lower frequencies. Furthermore, if it is sought to detect a class of 
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processes for which bandwidth ts proportional to the center frequency of the process, 


multiresolution techniques may provide the optimum detection technique. 


If a sequence s(n) 1s embedded in additive, white, Gaussian noise n(n), the resultant 
sequence 
x(n) = s(n) + N(n) (6.17) 
will also be a Gaussian random variable. If the noise process has n(n) a mean of zero then the 


mean of the signal plus noise x(n) 1s indicated by 


xia aes (en) (6.18) 


5. =o (6.19) 
In Chapter IV it was demonstrated that decomposition of a sequence by multiresolution 
structures constructed from cascaded perfect reconstruction quadrature muror filter banks is 
equivalent to projection of the sequence to be decomposed upon a space of orthonormal vectors. 


Decomposition by perfect reconstruction filter banks into subchannels, furthermore, can be 


interpreted as subdivision of the detail vector space W_, using the notation from Chapter IV, 


m? 
into another set of complementary set of subspaces es aa If, however, the filter banks do 


not satisfy the perfect reconstruction criteria, orthonormality of the equivalent projection 


Operation is not guaranteed. 


The the wavelet coefficients {b he -7 obtained from decomposing a sequence x(n) 


In, n 


were defined in Section IV.C as 
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bia = [x, Vma} (6.20) 
Now, if x(n) is a Gaussian random variable, then ban Will also be a Gaussian random variable. 


Furthermore, if the set of basis functions fy } is orthonormal, then the wavelet coefficients 


m,n 


ban Will be linearly independent. To obtain the expected value of b substitution of (6.17) 


m,n? 


into (6.20) produces 


be, Gye) = OE &. (6.21) 
Assuming the n has a mean of zero, then the mean of the wavelet coefficient b,, , for a signal 


s(n) embedded tn noise 1s simply the value of the wavelet coefficient evaluated in the absence of 


noise. 


To obtain the variance of b_ _, the squared mean 1s first 


m, n? 


Ciba =i Seana) sect Dl a (6.22) 


The expectation on the right-hand side of (6.22) is evaluated as 


EL(ns Vat = EX] nw) n(v) Wn. a(t) Ym.o(v) du dy}. (6.23) 
The expectation operator on the right-hand side of (6.23) can be passed through the double 


integral to operate on the terms n(u):N(v) producing 
E{n(u)-n(v)} =6,8(u - v) (6.24) 
where 6(t) denotes the Dirac delta function. By the integral sifting property of the Dirac delta 


function, (6.24) reduces to 


E{(n, Von. ee = 6, (Wan n? Vin, 3 a Ce (6.25) 
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The wavelet coefficients, therefore, simply comprise Gaussian random variables whose mean is 


; : ? 
(s, ae ) with a variance of Oe 


If detection of a known sequence in the presence of noise is required, one way to generate 
a decision statistic 1s to sum the squared magnitudes of all of the nodes in the transform space 
whose powers exceed a specified threshold. In performing multiresolution decomposition of a 


sequence, the essential information about a signal is contained in a linearly independent array 


D, ,, of non-uniform density such as that in figure 4.8. In detecting a known sequence in the 
presence of noise, it is conveniently possible to consider only the decomposition nodes in which 


a significant portion of the signal energy is represented. Mathematically, given a multiresolution 
decomposition represented by a lattice of points D, ,, the decision statistic can, equivalently, . 
constructed which excludes nodes contained in D; which are negligible in value. The subset of 
D; , of which the decision statistic is comprised can be limited to those nodes whose addresses 


are indicated by «: 
’ 2 2 
K={G,1): | rk | > p- max| Dik | 7 (6.26) 
J. 
Figures 6.2a, 6.2b, 6.3a and 6.3b present the results for statistics constructed from all 


nodes defined by (6.26) where p = 1/100. Table 6.2 indicates the mean square reconstruction 


error (3.32) which results from the value selected for 9. The information tn Table 6.2 ts 


intended as an indication of the consequences incurred from considering only the nodes in D, 
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where « is defined by (6.26). The decision statistics for the periodogram were generated in the 


same manner as for the multiresolution decompositions. 


TABLE 6.2--NORMALIZED MEAN-SQUARE RECONSTRUCTION ERROR RESULTING FROM 
TRUNCATED DECOMPOSITION. RECONSTRUCTION PERFORMED FROM SUBLATTICE DEFINED BY 
(6.26) WITH p = 1/100. EIGHT DECOMPOSITION STAGES WERE PERFORMED IN THE TWO-CHANNEL 
CASES WHILE FIVE-STAGE DECOMPOSITIONS ARE INDICATED FOR THE THREE-CHANNEL CASES. 


Test sequence (4.74a) Test Sequence (4.74b) | 


| Two-channel, 
| zero-subchannel 


| Two-channel, 
‘two-subchannel 


Two-channel, 
three-subchanne!l 


Two-channel, 
four-subchannel 


Three-channel, 
zero-subchannel 
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Figure 6.2a--Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
low-frequency transient (4.74a) embedded in noise with a -3 dB signal-noise ratio. Characteristics generated from 
500-realizations of 256-point random noise vector. 


For each technique indicated, 500 realizations of a 256-point, Gaussian random vector 


were generated. For each multiresolution technique reflected in Figures 6.2a and 6.2b, prior 
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generating the receiver operating characteristics, decision statistics were generated for separate 
noise realizations alone. The receiver operator characteristics for noise alone consisted of 
unity-slope straight lines. The decompositions were then performed for the test sequences 
(4.74a) and (4.74b) with additive noise and then for the noise vector alone. The decision 
statistics for the receiver operator characteristics consisted of all nodes in the composition whose 


magnitudes were greater than the one-percent threshold described above. 
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Figure 6.2b--Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
high-frequency transient (4.74b) embedded in noise with a -3 dB signal-noise ratio. Characteristics generated 
from 500-realizations of 256-point random noise vector. 


The receiver operating characteristics were, therefore, constructed from central and 
non-central chi-square distributions. The sum of the squares of the selected nodes of the noise 
decomposition, since they were generated from zero-mean noise, constitutes a chi-square 
distribution. Summing the selected nodes of the decomposition of signal plus noise produces a 


non-central chi-square distribution. 


In both Figures 6.2a and 6.2b, the spectrogram provided greater robustness than the 


multiresolution techniques. With the exception of the three-channel, zero-subchannel 
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decomposition applied to highpass test sequence (4.74b) in -3 dB signal-noise ratio, the 
differences were not, however, dramatic. Moreover, for test sequence (4.74b), the 
multiresolution curves were, with the named exception, closer to the performance of the 
spectrogram than for lowpass test sequence (4.74a). The greater disparity in the case of highpass 
test sequence (4.74b) may be attributable to the difference in resolution between the 


multiresolution and spectrogram techniques in the high-frequency regions of the spectrum. 


It is interesting to observe that, in Figure 6.2a, the best results among the multiresolution 
decompositions were obtained from the two-channel, zero-subchannel and two-channel, 
two-subchannel algorithms. Those two algorithms were constructed strictly from 
perfect-reconstruction filter banks and, therefore, produced the lowest reconstruction error. 
Additionally, in the case of highpass test sequence (4.74b), the best performance was obtained 
from the two-channel, three-subchannel decomposition. The two-channel, four-subchannel 


decomposition provided, by a small margin, the worst of all results. 
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Figure 6.3a--Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 


low-frequency transient (4.74b) embedded in noise with a -6 dB signal-noise ratio. Characteristics generated from 
500-realizations of 256-point random noise vector. 
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Figure 6.3b--Receiver operating characteristics for spectrogram and for various multiresolution decompositions of 
high-frequency transient (4.74a) embedded in noise with a -6 dB signal-noise ratio. Characteristics generated 
from 500-realizations of 256-point random noise vector. 


Results for a signal-noise ratio of -6 dB-are presented in Figures 6.3a and 6.3b. In the 
case of the lower signal-noise ratio, all of the multiresolution methods outperformed the 
spectrogram by a small measure. In the case of lowpass test sequence (4.74a) the best results 
were provided, as with the -3 dB signal-noise ratio test, by the two-channel, two-subchannel and 
the two-channel, zero-subchannel, respectively. This trend may be related to the superior 
reconstruction error provided by these techniques. In the case of highpass test signal (4.74b), on 
the other hand, the best results came from the four-subchannel and three-subchannel 
decompositions. This may have occurred as a result of the resolution of the decompositions. 
For the low-frequency transient, (4.74a), the three-subchannel and four-subchannel may have 
over-resolved the signal. For the high-frequency transient, (4.74b), the two-subchannel and 
zero-subchannel decompositions may have under-resolved the signal. Consequently, in the first 
case, the signal energy, for the higher resolution decompositions, may have been distributed over 


a greater-than-optimum number of time-scale bins. In the second case, the lower resolution 


Io 


decomposition may not have suffictently localized the stgnal to isolate tt from the effects of the 
noise. 


D. DECOMPOSITION OF STEADY-STATE HARMONICS AND FM CHIRPS 
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Figure 6.4a--Time-frequency plot of evolution of constituents of first 256-point test signal constructed from one 
steady-state harmonic and a quadratically chirped component. 
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Figure 6.46--Time-frequency plot of evolution of constituents of second 256-point test signal constructed from 
One steady-state hannonic and a quadratically chirped component. 


The time-scale resolutions of multiresolution decomposition techniques were somewhat 
addressed tn Section IV.E. However, in order to better understand the performance of these 


methods, a more controlled demonstration ts desirable. In thts section, results will be 
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demonstrated for a frequency-modulated, quadratic chirp converging to a steady-state harmonic 
component. The performance if the spectrogram method will first be illustrated. Then the 
results will be considered for two multiresolution techniques--the two-channel, four-subchannel 


and the two-channel, three-subchannel decompositions. 


Two test signals were generated for this demonstration. To illustrate the variations of 
resolution of scale provided by multiresolution decompositions in different regions of the 
frequency spectrum, it was desired to demonstrate the case in which the quadratic FM chirp 
converges from above to a steady-state harmonic in the lower half of the frequency spectrum and 
the case in which the FM chirp converges from below to a steady-state harmonic in the upper 
half of the frequency spectrum. Based on the experience of Chapter IV, it 1s expected that, for 


the first case, better resolution will be provided by the multiresolution techniques. 


The form for dynamic component of the test sequences is given by 
s(t) =cos(2-m- J, f(t) dt)+cos(2-n- fr- n) (6.27) 
where the instantaneous frequency f(t) possesses the form of a parabola whose vertex is, in the 
time-frequency plane, located at (t, f,), the coordinates of the final time and frequency in the test 
region. The time-frequency parabola satisfying this criterion assumes the form of 
f(t) = a-(t-t,)? + fp. (6.28) 
Evaluating the integral in the argument of the first cosine term on the right-hand side of (6.26) 


produces for the phase of the chirped component 


F(n) => - (inte)? -t? ) +f. (6.29) 
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To uniquely define f(t), it is necessary to specify the location of (t,, f,), the vertex of the 
parabola, and the location of some other point on the curve in the time-frequency plane. If the 
time t, = 0 and frequency f, of the chirp at the beginning of the test sequence are specified, the 


coefficient a from (6.28) becomes 


f.—f, 
a=. (6.30) 


Since 256-point test sequences are to be used, the final time tp = 256. It therefore remains to 


specify beginning and ending frequencies. 


For the first test sequence, a chirp which converges from below to a steady-state 
harmonic in the upper half of the frequency spectrum below the Nyquist frequency. The 


beginning and ending frequencies selected for the chirp were 


and 8 (6.3 Tap 


For the second test sequence in which the quadratic chirp converged from above to a steady-state 


harmonic in the lower half of the frequency spectrum, the frequencies selected were 


fp = <f, 


and (6.31b) 


The ideal, time-frequency evolution of the two test signals are plotted in Figures 4.6a and 4.6b, 


respectively. 


The spectrogram decompositions of the two test signals are plotted in Figures 6.5a and 


6.5b for two purposes. The first reason is to verify that the desired affects were obtained in 
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construction of the test sequences. The second reason for the periodograms is to provide a 
standard with which the performance of the multiresolution decompositions can be compared. 
The spectrogram plots were generated in a manner identical to that used in Section A of the 
present chapter. The data sequences were windowed by a 128-point Hanning window which 


was shifted in two-sample increments. Only the half of the frequency spectrum below the 


Nyquist frequency was plotted. 
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Figure 6.5a--Spectrogram decomposition of first 256-point test sequence constructed from a steady-state 
harmonic component and a frequency-modulated quadratic chirp. The decomposition was obtained by windowing 
the data sequence with a 128-point Hanning window shifted in increments of two samples. Only the half of the 
frequency spectruin below the Nyquist frequency is plotted. 


Examining Figures 6.5a and 6.5b, the two components of each test sequence can be seen 
to converge as expected. The steady-state harmonic produces a much better resolved indication 


than the chirped component. This occurs because of the smearing of the chirp which results 
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from the frequency shift occurring during the integration time of each FFT. Both 
decompositions represent the data with uniform resolution throughout the frequency spectrum. 


Consequently, the time, for both sequences, at which the components cannot be resolved is the 
same for both sequences at approximately the 230" sample. Consequently, from (6.28), the 


spectral resolution of the spectrogram decomposition is approximately 0.0043-f,. No attempt 
was made to optimize the window length with respect to resolution of the test signal 
components. The 128-point window provides good frequency resolution, however, because of 
the long integration time, the frequency of the chirped component has shifted significantly 


during each transform. 
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Figure 6.5b--Spectrogram decomposition of second 256-point test sequence constructed from a steady-state 
hannonic component and a frequency-modulated quadratic chirp. The decomposition was obtained by windowing 
the data sequence with a 128-point Hanning window shifted in increments of two samples. Only the half of the 
frequency spectrum below the Nyquist frequency Is plotted. 
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Figure 6, 6a--Two-channel, three-subchanne! multiresolution decomposition of first 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on (0, 13]. Subchannel decomposition was performed with 
pseudo-QMF filter bank designed in [11]. 


Examining Figures 6.5a and 6.5b, the two components of each test sequence can be seen 
to converge as expected. The steady-state harmonic produces a much better resolved indication 
than the chirped component. This occurs because of the smearing of the chirp which results 


from the frequency shift occurring during the integration time of each FFT. Both 
decompositions represent the data with uniform resolution throughout the frequency spectrum. 


Consequently, the time, for both sequences, at which the components cannot be resolved 1s the 
same for both sequences at approximately the 230" sample. Consequently, from (6.28), the 


spectral resolution of the spectrogram decomposition is approximately 0.0043-f.. No attempt 


was made to optimize the window length with respect to resolution of the test signal 


ba, 


components. The 128-point window provides good frequency resolution, however, because of 


the long integration time, the frequency of the chirped component has shifted significantly 


during each transform. 
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Figure 6.66--Two-channel, three-subchannel multiresolution decomposition of second 256-point test sequence 
constructed from steady-state harmonic component and quadratic clurp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was perfonned with 
pseudo-QMF filter bank designed in [11]. 


The two-channel, three-subchannel decomposition of the second test sequence 1s plotted 


in Figure 6.6b. In the case of this plot, the chirp becomes unresolvable from the steady-state 
harmonic at approximately the 230" sample. By (6.28), therefore, the frequency resolution is 


approximately 0.0043-f.. This number is quite close to the figure for the spectrogram 
decomposition. Furthermore, the two test sequence components become unresolvable around a 


scale of four which coincides detail subchannels obtained from the fourth stage of 
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decomposition. If the final frequency f, had occurred at a lower frequency, the convergence 


point between the two components would have occurred in the detail subchannels of a later 


decomposition stage resulting in even greater resolution. 
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Figure 6.7a--Two-channel, four-subchannel multiresolution decomposition of first 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on (0, 13]. Subchannel decomposition was perfonned with 
pseudo-QMF filter designed in Chapter V. 


Figures 6.6a and 6.6b present results of three-channel, four-subchannel multiresolution 
decomposition of the two test sequences. In Section IV.E, within the context of considering the 
decomposition of transients constructed from components of varying frequencies, the 
performance of the three-channel decomposition proved quite similar to that of the four-channel 
decomposition. This similarity appears in the case of decomposition of the present set of test 


sequences. In Figure 6.7a, the chirp and the harmonic become indistinguishable at 
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approximately sample number 108. This represents a slight improvement over the 
three-subchannel decomposition. In the case of Figure 6.7a, again using (6.28), the frequency 
resolution appears to be approximately 0.140-f,. This represents a subtle improvement over the 
performance of the three-channel decomposition. Nevertheless, in this region of the frequency 


spectrum, the spectrogram continues to provide performance which is superior by a significant 


margin. 
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Figure 6.7b--Two-channel, four-subchannel multiresolution decomposition of second 256-point test sequence 
constructed from steady-state harmonic component and quadratic chirp. Primary channel filters were based on 
Daubechies' orthonormal wavelet and scaling function on [0, 13]. Subchannel decomposition was performed with 
pseudo-QMF filter designed in Chapter V. 


Finally, the performance decomposition of the second test sequence is plotted in Figure 


6.7b. The chirp and the harmonic converge at approximately the same location. This result is 


also consistent the observations of Section IV.E. When decomposing transient test signals, it 
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was observed that increasing from three to four subchannels does not dramatically improve 
resolution of scale at large scales. The primary from the increase ts realized in the smaller 


scale--or higher frequency--regions. 
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VII. CONCLUSION 
A. SUMMARY OF RESULTS 


When compared with conventional periodogram techniques, the multiresolution techniques 
considered provided mixed results. In resolving signals at the extreme low end of the frequency 
spectrum below the Nyquist frequency, the multiresolution structures achieved their best results. 
In higher frequency regions, the multiresolution proved to be on par with or inferior to the 
periodogram. Overall, wavelet-like methods should be optimum for proportional bandwidth 
processes. The test signals employed performed marginally at demonstrating this strength in the 


upper regions of the frequency spectrum. 


With respect to other performance measures, the multiresolution techniques achieved, even when 
expanding fully by channel, decomposition of the test sequences with fewer calculations than the 
periodogram decomposition. In fairness to the latter method, the tmplementation demonstrated 
did not reflect any attempts to optimize the technique for the applications shown. However, in 
determining the number of calculations, the computational burden necessary to convert the 


resulting complex decomposition into a real representation was not considered either. 


The most interesting results were perhaps obtained from generation of the receiver operating 
characteristics. In the relatively high (-3 dB) signal-noise ratio demonstration the periodogram 
outperformed the multiresolution decomposition techniques by a small margin. However when 


the signal-noise ratio was decreased (-6 dB), the multiresolution decomposition achieved 
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superior results. Among the multiresolution structures, the results varied with the number of 


subchannels implemented and the region of the spectrum considered. 


Finally, in examining the application of multiresolution decomposition to dynamic signals, the 
results were entirely dependent on the region of the frequency spectrum in which the process 
being detected occurred. In the higher regions of the spectrum, the ability of the multiresolution 
techniques to resolve proximate narrowband components was poor at best. In the lower end of 
the spectrum, however, good results were achieved. 

B. RECOMMENDATIONS FOR FURTHER STUDY 

The most obvious direction in which to continue the study of multiresolution decomposition 
techniques 1s to improve the quality of the filter banks used for subchannel decomposition. In 
the ttre: implemented, decomposition into more than two subchannels was accomplished 
with pseudo-QMF banks. As indicated, methods for the design of perfect reconstruction filter 
banks of an arbitrary number of channels exist. The question arises as to whether results might 
not improve through the use of perfect reconstruction filter banks for decomposition into 


subchannels. 


Further study directed toward the determination of the optimum decomposition structure may 
also be warranted. In section 6.C, within the context of generation of comparing receiver 
operator characteristics for different structures, it was suggested that structures involving the 
decomposition of detail sequences into multiple subchannels may result in excessive resolution 
at low frequencies. As an alternative to decomposing the detail sequence into subchannels at 


each stage, limiting the performance of subchannel decomposition to the first two to three stages 
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may improve results. Alternatively, varying at each stage the number of subchannels into which 
the detail sequence 1s decomposed may also provide an option. For instance, a possible scheme 
would involve decomposition of the first stage detail sequence into three subchannels, the second 
stage detail sequence into two subchannels and performing no subchannel decomposition for the 


detail sequence at all subsequent stages. 


Additionally, the performance of multiresolution techniques for detecting dynamic signals in the 
presence of noise was not addressed. The periodogram decompositions of the quadratically 
chirped signals produced, through a spreading over multiple spectral bins, a significant reduction 
in the level of signal indication for the dynamic signal component. The steady-state harmonic 
component of the test sequence remained confined to one or two spectral bins while the dynamic 
component appeared spread over 12-16 spectral bins. Consequently, the signal levels of | 
equal-power dynamic and steady-state components of the signal varied by several decibels. This 
did not occur in the case of the multiresolution decomposition. Since, for the multiresolution 
decomposition the dynamic signal and the steady-state components remained equally well 
localized, the peaks of the dynamic signal component maintained the same approximate height 
as the steady-state component. It stands to reason, therefore, that the multiresolution detector 


may be less susceptible to noise when applied to the problem of a dynamic signal. 


Finally, further work is merited in the determination of the optimum manner in which to present 
multiresolution decompositions. For, a given signal, the essential information provided by 
multiresolution techniques is contained in a non-uniform lattice of nodes representing signal 


components spread, in varying degrees, over time. The method demonstrated in the present 
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work was to decompose the channel into channels and subchannels and then re-expand each 
channel individually in such a manner that the sum of the contents of each channel provides the 
reconstructed signal. This process of presentation introduces a significant computational burden. 
As was indicated in the section 6.B, the number of multiplications necessary to re-expand the 
signal always exceeded the number of computations for the actual decomposition by a 
significant margin. Brooks [33], on the other had presented only the coefficients of expansion in 
a surface plot with varying numbers of zeros inserted between nodes of information at different 
scales. This approach resulted in plots such as were presented in the present study for the a trous 
algorithm. Not only does this approach provide plots from which information 1s relatively 
difficult to extract, but as addressed in chapter 4, the issue of delays which are disparate among 
channels arises. Further consideration regarding the best manner in which to present the results 


of multiresolution decomposition 1s, therefore, merited. 
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APPENDIX A DETAILS OF DEVELOPMENT OF MALLAT'S 
ALGORITHM 


[n chapter 4, it was demonstrated that, since the scaling function vector 6 lies within the 


mtn 


span of the set of scaling function vectors, {, ,}, -z) any vector can be represented as a 


m+1,n 


Fourier series expansion in 9_, |: 


Om+i.n(t) = a (Ome Oa | Om.k(t). (A.1) 


Now the series coefficient, the inner product term contained in the summation (A.1), equals, by 


definition 
(Qmsi.n mk =| Om+i,n(t) " Om, k(t) dt. (m2) 


Substituting (4.14) into the mght-hand side of (A.2) produces 


Qremely2 2-m2 f (24m) .t—n)- (2 tk) di 
= —2™ f o(20™). t—n)- (2 - t-k) det | (A.3) 


v 


| O mei. n(t) . Dm, k(t) dt 





| rm | 


7m, 


2-™ f o(—S*) -9(2™ -t-k) dt 





ti 


j 
Now, (A.1) will become much simpler if the series coefficients can be expressed independent of 
the scale index m. Application of the change of variable of integration u = 2 "‘t - 2:n to the 


final line of (A.3) produces this result: 


Jadicern(t)eademme(t)idt =e J o(==**)-.g(2-™ -t-k) dt 
= > J o( = )-9((u+2-0)-k) du, (A.4) 


= > | o($)-9(u-(k-2-n) du 
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Since, therefore, by (4.14), Om y(t=27™? -O(2-" - t—n), (A.4) becomes 


= J o{$)-o(u-tk-2-n) | du = “aaa o.k-2n(u) du 





ex->) 
= [oro v5 Vo, 20 
Consequently, (4.19) 1s proved. From (A.5) and (4.19), (4.20) follows: 
Ometn(t) = = E (O10, Oo.kcon } Omlt) 
Developing (4.23) from (4.20) involves, first, substitution of the definition (4.14): 
2-im+1)2 O(2-lM*D . ty) = 2M (o: ibs 20 | o(2-™ -t—k), (A.6) 
Now, contained on each side of (A.6) are factors of 2°" and 2°"? Cancellation of those 
common factors produces 
o(2-) tn) ae (O10, bokre Oar aK). (A.7) 


Next, again introducing the change of variables u = 2° ""-t - 2:n results in 


o($] =2 (1.0, dorcza | o{u- (k- 2-n)). (A.8) 


Finally, translation of the index of summation k =k’ - 2:n in (A.8) yields 


o($)=E(o10, O04) Ou-W. (4.9) 


To generate the FIR filter defined in (4.23), it must first be specified that the scaling functions 


are compactly supported. Then, integrating both sides of (A.9) over their domains results in 


j o( +) au = E (O10 go) § uk) k) d (A.10) 
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Application of the changes in variables of integration v = w/2 and w =u - k to (A.10) produces 


om O(v) dv= - (or. On | J o(w) dw. (A.11) 


[f the common factors of | O(v) dv are removed from each side of (A.11), the result becomes 


2=5(O.0. Oink | (A.12) 


If h(k)=4 (O10. Oo , then 


Dante |: 
k 
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APPENDIX B MATLAB SOURCE CODE 


The source .\farlab code for primary algorithms employed during this study are provided in this 


appendix. 


A. LAPLACIAN PYRAMID 

% 

% This MATLAB function performs FULL-CHANNEL DECOMPOSITION 

% OP SEQUENGE IN ACCORDANCE WITH THE LAPMACIAN PYRAMID. 
% The Laplacian Pyramid algorithm 1s described in Michael 


% Unser, "An improved least squares Laplacian pyramid ... ", 
% Signal Processing 27 (1992) 187-203. 

% 

% The function syntax 1s: 

% P=lappmd(s, N) 

% where: "s" 1S a vector containing the sequence 

% to be decomposed. 

% "N" is an integer indicating the number 

% of levels of decomposition to be 

% performed. 

% "Pp" is an array whose rows contain the 

% sequence component corresponding 

% to that channel. : 

% The array "P" 1s computed such that summing along each 
2 column produces the reconstructed sequence. 


0 
Q 


function P=lappmd(s, N) 


Q 
O 


% The decomposition kernel is defined as "w." 

% 

w=[1/4-.2 1/4 .4 1/4 1/4-.2]; 

% 

% The input sequence "s" 1s reshaped as a row vector. 
% 

s=reshape(s,1,length(s)); 

% 

% The content 1s "s" is assigned to vector "pk1." The array 
% "D" ts initialized as a zero vector. 

% 

pkl=s; 

% 

% Decomposition of the sequence 1s performed. 


Wo 


0 


for k=1:N 
% 
% "pk" is reinitialized for next iteration. 
% 
pk=pk1: 
% The sequence 1s reduced using user-defined function REDUCE. 
pk 1=reduce(pk,w): 
% 
% "pk 1" is symmetrically windowed to one-half the length of "pk." 
% 
pkl=pk1(2:length(pk1)-1); 
% 
% The sequence "pk1" 1s expanded and the result assigned to "dk." 
% 
dk=expand(pk1,2*w); 
% | 
% "dk" 1s symmetrically windowed to the Iength of "pk." 
0% = 
dk=dk(3:length(dk)-2); 
% 


% The sequence "dk" is assigned to the kth row of array "D." The 
% length of "dk" 1s stored in the kth element of "L." 


% 
L(k)=min(length(dk),length(pk)); 
D(k,1:L(k))=pk(1:L(k))-dk(1:L(k)); 
end 
% 
% The final approximation ts assigned to the "k+1'th" row 
"Ay of oe 


O74 
70 


L(length(L)+1)=length(pk1); 
D(length(L),1:L(length(L)))=pk1; 


O7 

/0 

% The full-channel expansion ts calculated from the 

% rows of array "D" and the results assigned to the rows 
% of "P." The first row of "P" is equal to the first row 
Q/ eAN4 " 

/0 of D. 


Q/ 
40 


P( 1, =D): 
% 
for k=2:N+1 


192 


% The detail from channel "k" 1s assigned to vector "dk." 
% 
dk=D(k,1:L(k)); 
% 
°%) The channel content 1s expanded. 
a 
for j=2:k 
dk=expand(dk,2*w); 
dk=dk(3:length(dk)-2); 
end 
% 
% The expanded channel content ts assigned to row "k" of 
"%” array 'P." 
% 
P(k, 1:length(dk))=dk; 
end 


B. A TROUS ALGORITHM 

% 

% This MATLAB function calculates an approximation of the 
% DISCRETE WAVELET TRANSFORM of a function 

% using SHENSA'S ATROUS ALGORITHM discribed in "The 
% Discrete Wavelet Transform: Wedding the A Trous and 

“% Mallat Algorithms,” IEEE Transactions on Signal 

% Processing, Vol. 40, No. 10, Oct. 1992. 

% 

% The function syntax 1s: 

% W=shendwt(s, M, P_1, P_f) 

% where: "s" represents a vector of data to be transformed 


% "M" is an integer indicating the number of "voices" 
% tu be used to cover the frequency spectrum 

%) "W" is an array containing the coefficients of 

% the magnitude-squared wavelet transform 

%) Ch sy 

% "P_1" represents the first wavelet transform scale 
% to be calculated 

% "P_f" represents the final wavelet transform scale 
% The transform is calculated using a madulated Gaussian 
% window as an analyzing wavelet. If multiple voices are 


% specified, the projection of "s" on each voice will be 
% represented separately. The Lagrange, a trous interpolation 
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% filter used is obtained from convolving a four-point 
% Daubechies scaling filter with itself. 
function W=shendwt(s, M, P_i, P_f) 


(} 
0 


% First, the argument vector "s" 1s conditioned. If"s" 1s 
% defined as a column vector, it is converted to a row vector. 
% Secondly, "s" ts zeropadded to the next tnteger power of "2." 
% 
[rows,cols]=size(s); 
% 
if cols == | 
S=S.'; 
end 
clear rows. ciear Gols. 
% 
s=zeropad(s,1,2*ceil(log(length(s))/log(2))); 
% 


% Default values are imposed for starting and ending scales if 
% none are specified. 


% 

if exist(’P_1') == 
P_i=1; 

end 


if existe 
P FE floor(log(leneth(s))/log(2)): 
end 
% Ifthe number of voices ts not specified, a default value of 
% M=2 1s imposed. 


% 

if exist(‘M') == 
M=2; 

end 


Q 
Q 


% Next the analyzing wavelet must be calculated. The starting 

% point of this process 1s to specify the Gaussian window rolloff 
%) factor "beta" in accordance with the specified number of 

%) voices. (Shensa (6.31)). If M=1, the value of “beta” is defined 
% as "p1/(4*sqrt(2))." 


% 
if M == 1 
% 
beta=pi/(4*sqrt(2)); 
%) 
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else 


0 
0 


beta=1/(2*M): 

% 
end 
0 ) 
% The location of the center frequency "nu" ts assigned in accordance 
%) with Shensa (6.27). 
% 
nu=pi-sqrt(2) *beta; 
% 
% The voice scaling factor "a" 1s calculated using Shensa (6.32). 
% 
a=2’%(1/M); 

0 
% The region of support for the analyzing wavelet filter "g" is 
% approximated as the region for which the Gaussian window 
% is greater than 10%-3. Consequently, the filter impulse 
% response domain ts [-a*(M-1)*sqrt(14)/beta, a*(M-1)*sqrt(14)/beta]. 
% The length of "g" 1s represented by an odd integer. 
0) 0 
n=-ceil(a*(M-1)*sqrt(14)/beta):ceil(a*( M-1)*sqrt( 14)/beta); 
% 
% The analyzing wavelet is calculated for each voice in accordance 
% with Shensa (6.32). It is then normalized such that its peak 
% passband frequency response Is unity. 
% 
for k=1:M 
%o 

@(k,:)=exp(j *nu*(n/(a*(k-1)))). *exp(-(beta*2 *(n/a*(k-1)).%2)/2); 

% 
end 
eleal n; 
% 
% A variable "G" indicating the half-length of filter "g" is 
% defined. 
% 
G=ceil(length(g)/2); 
% 
% Next, the Lagrange interpolation filter is calculated. The filter 
% is obtained from "auto-convolution" of a Daubechie four-point 
“% DWT filter. 
% 
f=[1+sqrt(3) 3+sqrt(3) 3-sqrt(3) 1-sqrt(3)]/(4*sqrt(2)); 


ee 


% 

f=conv( f, fliplr(f)); 

f=f/sqrt(2); 

% 

% A variable "F" indicating the half-length of "f" 1s defined. 
%) 

F=ceil(length(f)2): 

% 

% The recursion 1s next executed according to Shensa (2.12a & b). 
Yo The sequence 1s filtered and decimated the first "P_i" times. 
% 


for k=1:P_1-1 
% 
s=conv(s,f); 
% 
s=s(F:2:length(s)-F); 
() 0 
end 
% 


% The output matrix "W" is initialized as a zero vector of 
% dimensions identical to those of "s." 

% 

W=zeros(s); 


0 
0 


for K—P Pee 
{) O 
% The data vector "s" 1s first filtered with each voice of "g.” 
% The squared magnitude of the result is assigned to the appropriate 
% column of "W." 
0 
0 
for n=(k-P_1)*M+1:(k-P_it+1)*M 
% 
% The row of "W" to be evaluated is initialized as zero. 
% 
W(n,:)=zeros( W(1,:)); 
% 
% The squared magnitude of the filter output is calculated 
% and assigned to "Wk." 
% 
Wk=sart(abs(conv(s,g(n-(k-P_1)*M,:))).%2); 
% 
% The elements of Wk" are assigned to the corresponding 
% columns of the row of "W" being calculated. 
% 
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W(n,2’(k-P_i):2%(k-P_1i):length( W))=Wk(G:length( Wk)-G+1); 


% 
end 
% 
% 
% "s" is filtered with the lagrange interpolation filter and then 
% @ecimated. 
% 
s=conv(s,f); 
% 
s=s(F:2:length(s)-F); 
% 
end 
GC. MULTIRESOLUTION DECOMPOSITIONS FROM CASCADES OF FILTER 
BANKS 
o 
% 
% This MATLAB function performs WAVELET-LIKE DECOMPOSITION 
% of a sequence using CASCADES OF QMF-TYPE FILTER BANKS 
% and then performs FULL EXPANSION BY CHANNEL 
% for each channel of decomposition. At each stage, the 
% sequence is decomposed into an approximation channel and 
% one or more detail channels. The detail channels can then 
% be further decomposed into subchannels. 
%) The function syntax 1s: 
% D=wavebank(s,h,N,g) 
% where: "s" is a vector representing the sequence to 
% be decomposed 
% "h" is an array whose rows contain the filter 
% coefficients for the primary analysis filter 
% bank in order of increasing center frequency. 
% "N" 1s an integer indicating the number of stages 
% of decomposition to perform. 
Ys g" is an array whose rows contain the filter 
“% coefficients for the analysis filter for the 
% subchannel decomposition in order of 
%) increasing center frequency. (OPTIONAL) 
% If "g" is not specified, the function, by default, does not 
% decompose the detail channel into subchannels. 
% 


oy 


function D=wavebank(s,h,N,g) 

% 

% The dimensions of "h" and "g" (if defined) are stored in 

% variables Re, Cc, Rs, Cs, indicating the array dimensions for 

Yo channel and subchannel decompositions, respectively. The number 


% of rows will be used as the decimation and expansion factors. 
% Then, the analysis filter banks are defined and stored in 

% “f for the primary channels and "1" for the subchannels. 

% Because of the effects of decimation, the synthesis filters 

“% must be scaled by a the decimation factor. 

% 


[Re,Cc]=size(h); 
f=Re*flipir(h); 
Le=length(h)-1; 


% 

Ls=0; 

Rs=1; 

if exist(‘g') == 1 

% 
[Rs, Cs]=size(g); 
1=Rs*fliplr(g); 
Ls=length(g)-1; 

% 

end “if 


() 
4) 


tt tf 


% The sequence "s" 1s reshaped as a row Vector. 

% 

s=reshape(s,1,length(s)); 

0%) 

% The sequence 1s decomposed. At each stage, the 
% approximation of "s" is stored in a vector "sk." 

% The next stage of decomposition is stored in "sk1." 


() 
0 


sk1(Re,:)=s; 


% 
for k=0:N-1 
o%/, 
% Sequence vector s Is reinitialized for next iteration. 
% 
sk=sk1; 
‘Yo 


% The sequence "s" 1s decomposed into primary channels by 
% the function REDUCE and the result Is assigned to "sk1." 


O 
0 
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for p=1:Re 


% 
dp=reduce(sk(Re,:),h(Re-ptl,:),Re): 
skl(p,1:length(dp))=dp; 

% 

end “op 

% 


“% "skl" 1s truncated to its proper post-"REDUCE" length. 


0 
M 


skl=sk1(:,l:length(dp)); 


% If "g" 1s defined, the sequence "dk!" is reduced into Re 
% subchannels. The lengths of the subchannel sequences for 
%o each channel are stored in a vector "L." 


for p=1:Rc-1 
if exist(‘g') == 
% 
for q=1:Rs 
dq=reduce(sk1(p,:),g(Rs-q+l,:),Rs); 
L(k+1)=length(dq); 
D((Re-1)*Rs*k+Rs*(p-1)+q, 1:L(k+1))=dq; 
% | 
end “oq 
% 
else 
% 
L(k+1)=length(dp); 
D((Re-1)*k+p,1:L(k+1))=sk1(p,:); 
% 
end “if 
end “op 
% 
end “%k 
% 
keyboard 
% 


% The sequence is reconstructed using full expansion by channel. 


QO 
0 


for k=1:N 

% 

% If "g" is defined, the detail "d" for each subchannel 1s 
% extracted from the corresponding row of "D." 
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{) 
MN 


if exist('g') == 


“A 


% 
% 
%, 
QO 0 


1 
0 


% 
Zo 
0 
0%, 


% 


% 


% 


Q 
Q 


QO 
0 


for p=1:Rc-! 
for q=1:Rs 


d=D(Rs*(Re-1)*(k-1)+Rs*(p-1)+q, 1: L(k)): 


The detail ts expanded with the corresponding subchannel 
synthesis filter. 


d=expand(d,i(Rs-qt+ l Ss 


The subchannel component ts expanded with the primary detail 
channel synthesis filter "f(@,:)" and then by the approximation 
channel synthesis filter "f(1,:)" the number of times 

appropriate for the number of stages of decomposition performed. 


d=expand(d,f(Rce+l-p,:),Rc); 
for m=1:k-1 

d=expand(d,f(1,:),Rc); 
end %“%m 


The total delay resulting from transmission through 
the system ts calculated and assigned to “delay.” It 
is then removed and the expanded sequence "d" ts 

reintroduced to the corresponding row of array "D." 


delay=sum((Rc).*(0:k-1])*Le+Re“*k*Ls; 
d=d(delayt+|:delay+length(s)); 
M(k)=length(d); 
D(Rs*(Re-1)*(k-1)+Rs*(p-1)+q, 1:M(k))=d; 


end %q 
end “op 


If "g" is not defined, the subchannel reconstruction steps 


are bypassed. Reconstruction of the primary detail channel 
is performed. 
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[Grp lec | 

d=D((Re-1)*(k-1)+p, |:L(k)): 
“% 
% The detail sequence is expanded by the primary detail channel 
% synthesis filter "f(Rc+l-p,:)" and then by the approximation 


% channel synthesis filter "f(1,:)" as many times as 
% appropriate for the number of stages of decomposition performed. 
% 

d=expand(d,f(Rc+1-p,:),Rc); 
% 

for m=1:k-1 
% 

d=expand(d,f(1,:),Rc); 

“% 

end %m 
% 


% The total delay resulting from transmission through 
% the system ts calculated and assigned to "delay." It 
% is then removed and the expanded sequence "d" is 
% reintroduced to the corresponding row of array "D." 
% 
delay=sum((Rc).*[0:k-1])*Le; 
d=d(delay+1:delay+length(s)); 
M(k)=length(d); 
D((Re-1)*(k-1)+p, 1:M(k))=d; 
% 
end “%p 
% 
end “if 


( 
My 


end %k 

“% 

‘ The final approximation sequence Is selected as the final 
% row of "sk1." 

% 

skl=sk1(Re,:); 

% 

% The final approximation signal is expanded with primary 
% approximation channel synthesis filter according to the 
% number of decomposition stages. 

%) 

for m=1:N 


0 
) 
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sk l=expand(sk1,f(1,:),Re): 
M%) 
end %m 
% 
% The expanded approximation sequence is added as the final row 
“% of alta y wie 
% 
“% The total delay resulting from transmission through 
% the system is calculated and assigned to “delay.” It 
% 1s then compensated and the expanded sequence "sk1” is 
% reintroduced to the final row of array "D." 


Q 
VU 


delay=sum((Rc).*[0:N-1])*Le; 

sk l=sk1(delay+1:delay+length(s)); 
M(length(M)+1)=length(sk1); 
D(Rs*(Re-1)*N+1,1:M(length(M)))=sk1; 


D. PERIODOGRAM DECOMPOSITION 


% This MATLAB function calculates the PERIODOGRAM SPECTRAL ESTIMATE 
% for a specified input data vector. The function syntax 1s 


% 

% S=pgram(s,w, N); 

% 

% where: "s" 1s the data vector to be analyzed. 

% "w" represents a window with which data will be windowed 
% "N" is an integer indicating the number of steps to 

% shift the window for each step 

% (default = length(w)/2) 


% "S" 1s a perlodogram-type output matrix. 

{) 0 

% The argument data vector 1s subdivided with 50% overlap. Each 
% row of the output array represents the transform of a subdivision 
% of the original data vector. 

Q 0 

function S=pgram(s,w,N); 

% 

% Check to see if "w" is defined. If '"w" is not defined, 

% a rectangular window of length "length(s)” 1s imposed as 

% a default. 

“A 

if exist(‘w') == 
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% 


() 
0 


end 
vy 


() 
0 


”) 


W=ones(s); 


A default vale is imposed for "N" 


if exist("N') == 0 


%) 
% 
end 
% 
ie) 0 


() 
Q) 


Yo 


N=length(w)/2; 


Convert column vector data vector arguments to row vector 
format. 


s=reshape(s,length(s),1); 
w=reshape( w,length(w), 1); 


% 
% 
Yo 


"s" 1s zeropadded to an integer number of window lengths. 


for k=1:(length(s)-length(w))/N+1 


% 


“) 


Q 
0 


end 


0 
0 


) 
0 


% 


{ 
”) 


r=fft(w. *s((k-1)*N+1:(k-1)*N+length(w))); 


S(.,k)=r(1:length(w)/2); 


In accordance with the definition of the periodogram, the 
periodogram magnitude 1s divided by the window length. 


S=S/length(w); 


EB. 
% 
) 0 
% 
% 
() 0 
() O 
EG 
“% 
% 


This function expands a two-scale, difference equation of 
the form 
phi(x)=sum(c_k*phi(2*x-k)). 
The function syntax is 
[phi, psi]=twoscale(c,P,eigtol) 
where the function arguements are: 
"cis a coefficient vector [c_O c_l... c_N-1] 
"P" 1s the number of dyadic points evaluated between 
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% each integer in the domain. 


% “eigtol" is indicates the maximum deviation from unity 
% which will be tolerated in calculating the eigenvalue 
% corresponding to the eigenvector containing the 

% values of "phi" for integer arguments. 

% (DEFAULT = 1E-3) 

“% The output variables are 

% "phi" represents the scaling function satisfying the 

"% difference equation 

% "psi" 1s the corresponding wavelet. 

% 

function [phi,psi]=twoscale(c,P, eigtol) 

% 

% If "eigtol” is not defined, a default vale is imposed. 

% 


if exist(‘eigtol’) ~= | 
eigtol=le-3; 


end 

% The argument vector "c" is reshaped as a column vector. 

% 7 
c=reshape(c,length(c),1); 

% 

% The condition that the sum of the elements of "c" equals "2" 
% 1s Imposed. 

Yo 

c=2* c/sum(c); 

% 

% First the scale function is evaluated for integer points 

%, on its domain. A zero-padded version "c_n" of the coefficent 
% vector "c" 1s created. 

% 

c_n=[zeros(l:length(c)-1)'; c; zeros(1:length(c)-1)']; 

0 0 

% A matrix "C" containing values of the coefficient vector 

%) "c" 1s generated. 


% 
for k=1:length(c) 


0 
Q 


C(k,:)=flipud(c_n(2*k-1:2*(k-1)+length(c)))’; 


0 
0 


end 

% 

% The matrix "C" is truncated to eliminate rows corresponding 
% to the trivial solutions corresponding to non-unity 
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% wen and C_N-|”. 


% 
if c(1) ~= 1 
% 
C=C(2:length(C),2:length(C)); 
07) 
end 
10 
if c(length(c)) ~= 1 
%) 
C=C(1:length(C)-1,1l:length(C)-1); 
% 
end 
% 
% The resulting matrix 1s checked to ensure that it contains 
% at least one eigenvalue of unity. If one unity eigenvalue 
% 1S present, the function continues. Otherwise, it terminates 
% In an error message. 
% 
% The eigenvectors and eigenvalues of "C" are calculated. 
% The diagonal matrix containing the eigenvalues is converted to 
% a column vector of eigenvalues. 
% 
ey |=cig(C); 
D=diag(D); 
% 


if min(abs(D-ones(D))) <= eigtol 


0 
0) 


%) The address of the unit eigenvalue is determine. The 
% elgenector corresponding to that address is defined as "phi." 
%) Values of "zero" are imposed on the endpoints for non-unity 
%) "Cal OL, coN-1", 
% 
K=find(abs(D-ones(D)) == min(abs(D-ones( D)))); 
R(T; 
phi=E(:,K); 
vi, 
ie t= | 
% 


phi=[0;phi]; 
end 
% 
if c(length(c)) ~= 1 


” 
end 
% 
“% 
% 
%) The values of the scaling function are caluclated recursively 
“) for dyadic values of "phi." The values of “phi” at half-integer 
% domain points 1s merely the convolution of "phi" and "c." 
% 
N=length(c)-1; 
for k=0:max(0,P-1) 
% 
% A circulant, convolution matrix "M" is created. 
% 
% The matrix "M" is cleared from memory for the next iteration. 
% 
clear M 
% 
for m=0:N 
% 
M(:,m+1)=[zeros(m*2*k, |);phi;zeros((N-m)*2*k, 1)]; 
Q 0 
end 
% 
%size(M) 
% The convolution ts evaluated. 
% 
phi=M *c; 
% 
% 
end 
% 
% The corresponding wavelet 1s calculated. 
% 
n=0:length(c)-1; 
c=flipud(c).*((-1).4n’); 
psi=M*c; 
% 
% An error message 1s displayed in the event no 
% eigenvalues of unity are present. 
% 
else 


phi=[ph1,;0}; 
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% 

‘Filter string produces no unity-value eigenvalues.’ 
Y% 
end 


Ae GENERATION OF RECEIVER OPERATING CHARACTERISTICS 


% This AfA4 TLAB script file generates receiver operating characteristics 
% for various multiresolution decompositions, saves the results 

% to files, and when complete, terminates \/47ZLAB. 

% 

“% Generate roc curve for two-channel, four-subchannel, 

% SNR = - 3 dB, high-frequency transient 

% 

load sIf 

load pqmfb4 

g4=h; 


load wavebank 

D=wbhdec(s,h,8,24); 

K=find(D.*2 >= max(max(D.*%2))/100); 
rand(‘normal’); 

sig=sqrt(mean(s.*2)); 


% 

‘two-channel, four-subchannel’ 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8,¢4); 
HO0(k)=sum(sum(D(K).%2)); 
D=wbdec(nts,h,8,24); 
H1(k)=sum(sum(D(K).%2)); 

end 

% 


(PFA231f,PD231f]=roc( HO, H1, 50); 

save roc23lf PFA23lf PD23lf 

% 

clear 

0% 

clear 

% 

% Generate roc curve for two-channel, three-subchannel, 
% SNR = - 3 dB, high-frequency transient 

% 
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ie adlusiit: 

load wavebank 

D=wbdec(s,h,8.g); 

K=find(D.*2 >= max(max(D.%2))/100); 

rand(‘normal’): 

sig=sqrt(mean(s.*2)): 

% 

‘two-channel, three-subchannel' 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8,g); 
H0O(k)=sum(sum( D(K).*2)); 
D=wbdec(n+s,h,8,g); 
H1(k)=sum(sum(D(K).%2)); 

end 

% 

(PFA231f,PD23ln|=roctnGmiat, Su), 

save roc23lf PFA23If PD23lf 

% 

clear 

%) 

% Generate roc curve for two-channel, two-subchannel, 

% SNR = - 3 dB, high-frequency transient 

% 

load slf 

load wavebank 

D=wbdec(s,h,8,h); 

K=find(D.*%2 >= max(max(D.%2))/100); 

rand(‘normal'); 

sig=sqrt(mean(s.*2)); 


%) 

‘two-channel, two-subchannel’ 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8,h); 
HO(k)=sum(sum(D(K).%2)); 
D=wbdec(nts,h,8,h); 
H1(k)=sum(sum(D(K).%2)); 

end 

% 


[PFA22iePB22\f]=roc(HOerns 56y, 
save roc22If PFA22If PD22If 


) 
0 
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clear 
% 
4) 
0 
a Generate roc curve for two-channel, zero-subchannel, 
%) SNR = - 3 dB, high-frequency transient 


0 
My 


load sif 

load wavebank 

D=wbdec(s,h,8); 

K=find(D.*2 >= max(max(D.%2))/100); 

rand(‘normal'); 

sig=sqrt(mean(s.*2)); 

‘two-channel, zero-subchannel' 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8); 
HO0(k)=sum(sum(D(K).%2)); 

=wbdec(n+s,h,8); 

H1(k)=sum(sum(D(K).‘2)); 

end 

% 

(PFA20If,PD201f]=roc( HO, H1, 50); 

save roc20If PFA20If PD20If 


(0) 
i) 


clear 


0 
0 


% 
%) Generate roc curve for three-channel, zero-subchannel, 
% SNR = - 3 dB, high-frequency transient 
% 
load slf 
load wavebank 
D=wbdec(s,g,8); 
K=find(D.*2 >= max(max(D.%*2))/100); 
rand(‘normal’); 
sig=sqrt(mean(s.*%2)); 
% 
‘three-channel, zero-subchannel' 
for k=1:500 
n=4*sig*rand(s); 
D=wbhdec(n,g,8); 
H0(k)=sum(sum(D(K).%2)); 
D=wbdec(n+s,g,8); 
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H1(k)=sum(sum(D(K).%2)); 
end 
%) 
(PFA30If, PD30if]=roc( HC, 1 5G}, 
save roc30If PFA30If PD30If 
% 
clear 
clear 
% 
% Generate roc curve for two-channel, four-subchannel, 
% SNR = - 3 dB, high-frequency transient 
%) 
load shf 
load pqmfb4 
at 
load wavebank 
D=wbhdec(s,h,8,24); 
K=find(D.*2 >= max(max(D.%2))/100); 
rand(‘normal'); 
sig=sqrt(mean(s.*%2)); 


% 

‘two-channel, four-subchannel'’ 

for k=1:500 © 
n=4*sig*rand(s); 
D=wbdec(n,h,8,2¢4); 
HO0(k)=sum(sum(D(K).%2)); 
D=wbhdec(n+s,h,8,¢4); 
H1(k)=sum(sum(D(K).%2)); 

end 

% 


([PFA23hf,PD23hf]=roc(HO, H1, 50); 

save roc23hf PFA23hf PD23hf 

% 

clear 

%) Generate roc curve for two-channel, two-subchannel, 
% SNR = - 3 dB, high-frequency transient 
% 

load shf 

load wavebank 

D=wbdec(s,h,8,h); 

K=find(D.*%2 >= max(max(D.%2))/100); 
rand(‘normal'); 
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sig=sqrt(mean(s.*2)); 

% 

‘two-channel, two-subchannel’ 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8,h); 
HO(k)=sum(sum(D(K).%2)); 
D=wbdec(n+s,h,8,h); 
H1(k)=sum(sum(D(K).%2)); 

end 

% 

[PFA22hf,PD22hf]=roc(HO, H1, 50); 

save roc22hf PFA22hf PD22hf 

% 

clear 

% 

% Generate roc curve for two-channel, zero-subchannel, 

% SNR = - 3 dB, high-frequency transient 

% 

load shf 

load wavebank 

D=wbdec(s,h,8); 

K=find(D.*2 >= max(max(D.%2))/100); 

rand(‘normal’); 

sig=sqrt(mean(s.*2)); 


% 

‘two-channel, zero-subchannel’ 

for k=1:500 
n=4*sig*rand(s); 
D=wbdec(n,h,8); 
HO(k)=sum(sum(D(K).%2)); 
D=wbdec(nts,h,8); 
H 1(k)=sum(sum(D(K).%2)); 

end 

% 


[PFA20hf,PD20hf]=roc(HO, H1, 50); 
save roc20hf PFA20hf PD20hf 

% 

%) 

quit 


0 
0 


‘ce GENERAL-PURPOSE ROUTINES CALLED BY OTHER ROUTINES 


% This function zero-pads an argument matrix forming 

% a new matrix of the specified size. The original matrix 

% will be in the upper left-hand corner the new matrix. 

% 

% The function syntax 1s: 

% 

% X=zeropad(x,N1,N2) 

% 

% where: 

% "x" 1s the initial argument matrix 

% "NI" 1s the total number of rows for the new matrix 
% "N2" is the total number of columns for the new matrix 
% 

function X=zeropad(x,N1,N2) 

% 

[M1,M2]=size(x); 

% 

X=[x zeros(M1,N2-M2)]; 

% 


X=[X; zeros(N1-M1,N2)]; 


% 

% This MATLAB function reduceS A SEQUENCE IN THE SENSE OF 

% THE EXPAND OPERATION EMPLOYED IN LAPLACIAN PYRAMID 
% DECOMPOSITION. The operation entails FIR filtering 


%) followed by decimation in time. 

% 

% The function syntax 1s 

% fk 1=reduce(fk, h, M) 

% where: fk is a vector containing the sequence to be expanded 
% h is a vector containing the FIR filter coefficients 
“% M is an integer indicating the decimation 

% factor (default == 2) 

a 

% Regardless of the shape of the input vectors, output vectors 

% are returned as row vectors. 

% 

function fkl=reduce(fk, h, M); 

% 

% A default value of two 1s imposed on "M." 

% 


to 
— 
tr 


If exist('M') ~= | 


M=2: 
end 
% 
% The sequences "fk" and "h" are reshaped as row vectors. 


0 

fk=reshape( fk, 1 length(fk)); 

h=reshape(h,1,length(h)); 

() 0 

% Mnesequence ik 1s filteredwith hand the result is 
% assigned to "fkl." 

% 

fkl=conv(fk,h); 

% 

% The sequence "fk1" 1s decimated by a factor of "M." 
% 

fkl=fk1(1:M:length(fk1)); 


% 

% This MATLAB function EXPANDS A SEQUENCE IN THE SENSE OF 
% THE EXPAND OPERATION EMPLOYED IN LAPLACIAN PYRAMID 
% DECOMPOSITION. The operation entails FIR filtering 


% followed by decimation in time. 
% | 
% The function syntax is 
% fk 1=expand(fk, h, M) 
% where: fk 1S a vector containing the sequence to be expanded 
% h 1s a vector containing the FIR filter coefficients 
% M is an integer indicating the decimation 
% factor (default == 2) 
% 
% Regardless of the shape of the input vectors, output vectors 
% are returned as row vectors. 
% 
function fkl=expand(fk, h, M); 
% 
% A default value of two is imposed on "M." 
% 
if exist('M’) ~= | 
M=2; 
end 
% 
% The sequences "fk" and "h" are reshaped as row vectors. 
% 
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fk=reshape( fk, | ,length(fk)); 
h=reshape(h, | ,length(h)); 


% 

% Zero-insertion 1s performed on the input sequence "fk." 
% The result 1s assigned to "fk." 

% 


fk l=zeros(1,M{*length(fk)-1); 

tk1(1:M:length(fk1))=fk; 

% The sequence "fk" 1s filtered with "h" and the result is 
% assigned to "fk1." 

% 


fk 1l=conv(fk1,h); 


% 

% This function performs numerical integration of a 

% vector in accordance with Simpson's Rule. The function 
% arguments are: 

% 

% "x"--the domain of the argument function 

% "y"--the range of the argument function. 

% 

% The function's syntax 1s: 

% 

% Y=simpson(x,y) 

% 

% A regular partition 1s assumed. 

% 

function Y=simpson(x,y) 

% 

% Vectors are checked for row- or column-vector format. Column 
% vectors are converted to row vectors. 

“% 


dimens=size(x); 
N=max(dimens); 
if dimens(1) ==N, 


X=x'5 
y=’; 
end 
% 
% The partition width is determined and the argument range vector 
% is multiplied by the partition width. 
% 


y=y *((max(x)-min(x))/max(size(x))); 


214 


% The integral 1s evaluated. 

% 

mi) —y(1); 

for k=2:N 
Y(k)=v(1:k)*ones(y(1:k)’); 

end 

%) 

% 


% Given two, equal-length vectors, this MATLAB function 
% GENERATES HISTOGRAMS AND CONVERTS THEM TO RECEIVER 
% Oren A TING CHARACTERISTICS. 


% 

% The function syntax 1s: 

% [PFA, PD]=roc(H0, H1) 

% where: "HO" is a vector containing "hypothesis false” 
% realizations of a random process. 

% "HI" is a vector containing "hypothesis true” 
% realizations of a random process. 

% "PFA" tsa vector representing the probability 
% of false alarm 

% "PD" is a vector representing the probability 
% of detection. 


) 0 
function [PFA, PD]=roc(HO, H!, Nobins) 
% 
% The minimum and maximum bin locations are identified. 
% 
Bmin=min(min({(H0;H1])); 
Bmax=max(max({H0;H1])); 
% 
% A vector 1s created representing the bin locations. 
Va 
if exist("Nobins’) == 
% 
Nobins=length(H0)/25; 
% 
end 
B=Bmin:(Bmax-Bmin)/Nobins:Bmax; 
% 
% The MATLAB function "HIST" ts invoked to produce histograms 
% of the realization vectors. 
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% 

[pfa,X]=hist(HO,B); 

[pd,X]=hist(H1,B); 

07, 

‘“% The histograms are integrated using function "SIMPSON." 
PF A=simpson( X,pfa); 

PD=simpson( X,pd); 

% 

% The integrals are normalized for maximum values of unity. 
% 

PFA=PFA/max(PFA); 

PD=PD/max( PD); 

% 

% The actual probability detect/probability false-alarm 

% plots are the complements of the probability distribution 
% functions. 


Yo 
PFA=1-PFA; 
Abe) 
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